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We implemented the charge self-consistent combination of Density Functional Theory and Dy- 
namical Mean Field Theory (DMFT) in two full-potential methods, the Augmented Plane Wave 
and the Linear Muffin- Tin Orbital methods. We categorize the commonly used projection methods 
in terms of the causality of the resulting DMFT equations and the amount of partial spectral weight 
retained. The detailed flow of the Dynamical Mean Field algorithm is described, including the 
computation of response functions such as transport coefficients. We discuss the implementation of 
the impurity solvers based on hybridization expansion and an analytic continuation method for self- 
energy. We also derive the formalism for the bold continuous time quantum Monte Carlo method. 
We test our method on a classic problem in strongly correlated physics, the isostructural transition 
in Ce metal. We apply our method to the class of heavy fermion materials Celrlns, CeCoIns and 
CeRhlns and show that the Ce 4/ electrons are more localized in CeRhlns than in the other two, 
a result corroborated by experiment. We show that Celrlns is the most itinerant and has a very 
anisotropic hybridization, pointing mostly towards the out-of-plane In atoms. In CeRhlns we sta- 
bilized the antiferromagnetic DMFT solution below 3 K, in close agreement with the experimental 
Neel temperature. 

PACS numbers: 71.27.+a,71.30.+h 



I. INTRODUCTION 

One of the most active areas of condensed matter the- 
ory is the development of new algorithms to simulate and 
predict the behavior of materials exhibiting strong cor- 
relations. Recent developments in the dynamical mean- 
field theory (DMFT^, a powerful many-body approach, 
hold great promise for more accurate and realistic de- 
scriptions of physical properties of this challenging class 
of materials. 

The crucial step towards realistic description of 
strongly correlated materials was the formulation of 
DFT+DMFT^El, a method formed by the combination 
of density functional theory (DFT) and DMFT (for a 
review see Ref. [5]). To date, this method already has 
substantially advanced our understanding of the physics 
of the Mott transition in real materials and demonstrated 
its ability to explain phenomena including the structural 
phase diagr ams o f actinides^Sl, phonon respons e^, opt ical 
conductivity^^, valence and x-ray absorptiorPHSl &n d 
transport!^ of archetypal strongly correlated materials. 

At present, much effort is devoted to the development 
of a robust and precise implementation of DFT+DMFT 
using state of the art DFT electronic structure codes^" 2 ^ 
and advanced impurity solveriPfflHI This article de- 
scribes in detail the implementation of this method 
within full-potential codes. There are three major is- 
sues that arise in DFT+DMFT implementations: i) qual- 
ity of the basis set, ii) quality of the impurity solvers, 
and iii) choice of correlated orbitals onto which the full 
Green's function is projected. Modern DFT implemen- 
tations largely resolve the first issue, recent development 
of new impurity solvers^^ 12 ^ 127 1 29 ' have focused atten- 
tion on the second, while the third is rarely discussed in 



the literature. Many DFT+DMFT proposals in the lit- 
erature are based on downfolding to low energy model 
Hamiltonian d 2 ! 19 ! 20 ! 28 !, which requires an atomic set of 
orbitals and treats the kinetic operator on the level of 
an effective tight binding model. In contrast, we avoid 
the ambiguities of downfolding and instead keep the ki- 
netic part of the Hamiltonian and electronic charge ex- 
pressed in a highly accurate full potential basis set. The 
advantage of our method is its ability to perform fully 
self-consistent electronic charge calculations. We con- 
centrate here on the Linear Augmented Plane Wave basis 
(LAPW^as implemented in the Wien2K code^and the 
LMTO basis as implemented in LmtArlP2l, in combina- 
tion with the imp urity solvers based on the hybridization 
expansio n 1 21 1 23 ^ . 

The first half of the article introduces the basic steps of 
implementing the DFT+DMFT algorithm and provides 
a pedagogical introduction to the method. Section [IT] is 
devoted to a crucial element of the DFT+DMFT formal- 
ism, namely the projection of the full electronic Green's 
function to the correlated subset. We show that the pro- 
jection used in the LDA+U method leads to non-causal 
DFT+DMFT equations, while the projection on to the 
solution of the Schrodinger equation within the Mufhn- 
Tin (MT) spheres misses electronic spectral weight. We 
propose a new projection that leads to causal DMFT 
equations and captures all electronic spectral weight. 
Section III derives the DFT+DMFT equations from a 
Baym-Kadanoff-like functional formalism. Section |IV| 
provides a detailed flowchart of all the steps of the al- 
gorithm. In section [V] we discuss the necessary changes 
to the tetrahedron method when used in the context of 
DMFT. Section |VI| described the algorithm to compute 
transport properties within DFT+DMFT. Section |VII| 
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describes the impurity solvers based on the hybridiza- 
tion expansion, the One Crossing Approximation (OCA), 
and the Bold continuous time quantum Monte Carlo al- 
gorithm (6-CTQMC). Finally, section VIII discusses a 
new algorithm for analytic continuation of the self-energy 
from the imaginary to real axis. 

In the second half of the article, we describe the re- 
sults obtained by applying our new implementation of 
DFT+DMFT to several correlated materials. As a first 
test of the algorithm, in section |IX| we present its appli- 
cation to elemental cerium. Section [X] is devoted to a 
class of heavy fermion materials, CeRhlns, CeCoIns and 
Celrln 5 , dubbed Ce-155 materials. We show the differ- 
ence in the electronic structure among these three mate- 
rials and demonstrate that the Ce 4/ electrons are most 
localized in CeRhlns and order antiferromagnetically be- 
low T/v ~ 3K, in agreement with experiment, while the 
Ce 4/ electrons are most itinerant in Celrln5. We ex- 
plain the origin of the subtle difference between the three 
Ce-115 compounds from the electronic structure point of 
view. 



II. PROJECTION ON TO CORRELATED 
ORBITALS WITHIN FULL-POTENTIAL 
METHODS 

DFT+DMFT contains some aspects of band theory, 
adding a "frequency-dependent local potential" to the 
Kohn-Sham Hamiltonian. It also contains some aspects 
of quantum chemistry, carrying out an exact local con- 
figuration interaction procedure by summing all local di- 
agrams, which requires the definition of an "atomic-like" 
or "local" Green's function. The operation of extracting 
the local Green's function Q(r,r') from the full Green's 
function G(r, r') is called projection (or truncation). The 
reverse operation of expressing the local time-dependent 
potential derived from the solution of the atomic 

problem in the presence of a mean- field environment, is 
called embedding. The various DFT+DMFT implemen- 
tations differ not only in the choice of basis set, but also 
in the choice of the projection-embedding step. These 
ingredients are sketched schematically in Fig. [T] The 
projection-embedding step connects the atomic and solid 
state physics, and its proper definition is a conceptual is- 
sue of DFT+DMFT method. 

In the current formulation of DFT+DMFT^ 5 ' 33 ' 34 !, 
one must define the correlated orbitals to which 
the Coulomb correlation is applied, i.e., £(r, r') = 
Ylff> Xf( r )^'??'X|'( r ')) where x?( r ) i s a localized or- 
bital. Usually, this is achieved by transforming 
the DFT Hamiltonian to a set of localized Wannier 
orbitals. These Wannier orbitals are then identi- 
fied as the local correlated orbitals of DMFT. Var- 
ious choices of these orbitals were propos ed in the 
literature, including tight-binding LMTO'a^, non- 
orthogonal LMTO'sP, Nth-order Muffin-Tin orbitalfP, 
numerically-orthogonalized LMTO's 3 ^, and maximally- 
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FIG. 1: Schematic diagram of the projection-embedding 
step in the DFT+DMFT algorithm. The full Green's func- 
tion of the solid G(r, r') is truncated to its local counterpart 
PG — Gll'- The impurity solution delivers an effective local 
potential, which is embedded (E) into the Dyson equation 
of the solid. The DMFT self-consistency condition (DMFT- 
SCC) connects the two. 



localized Wannier orbitals^EH. The basis functions must 
fully respect the symmetries of the problem and be atom- 
centered, rather than bond-centered. Hence maximally- 
localized Wannier functions^ are not a good starting 
point for DMFT. 

Localized basis sets are a better starting point for our 
purposes, but the non-orthogonality of these sets pose 
a serious challenge. Straighforward orthogonalization 
mixes the character of the orbitals, resulting in mixed the 
partial occupancies and partial density of states, leading 
to incorrect partial electron counts. For example, within 
modern DFT implementations, cerium metal has approx- 
imately one 4/ electron. Naive orthogonalization results 
in a considerably higher 4/ electron count, leading to an 
unphysical DMFT solution. 

Even more challenging is the formulation of the good 
localized orbitals in full-potential basis sets. Here, mul- 
tiple basis functions are used to obtain more variational 
freedom. To implement DMFT in such basis sets, the 
group of orbitals representing the correlated electrons 
in the solid must be contracted to form a single set of 
atomic-like heavy orbitals, i.e., one 4/ orbital per Ce 
atom, one 3d orbital per Fe atom, etc. 

A straighforward projection on to the orbital angular 
momentum eigenfunctions Yi m = Yl leads to non-causal 
DMFT equations, which result in an unphysical auxiliary 
impurity problem. The second often-employed choice 
is the projection on to the solution of the Schrddinger 
(Dirac) equation inside the MT sphere Ui{E u , t)Yl{v). 
While this choice is certainly superior to the straighfor- 
ward projection, it does not take into account the contri- 
butions due to the energy derivative of the radial wave 
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function ui{E v , r)Y^(f) and the localized orbitals (LO) 
at other energies ui{E' v ,r)Yi,{r), and hence misses some 
electronic spectral weight of the correlated orbital. Alter- 
native choices are possible which simultaneously capture 
all spectral weight and obey causality. We implemented 
one of them and we believe it is superior to other choices 
in the literature. 

The central objects of DMFT are the local Green's 
function and the local self-energy of the orbitals within 
the correlated subset. We specify the projection scheme 
by the projection operator P(rr' ,tLL'), which defines 
the mapping between real-space objects and their orbital 
counterparts (r, r') — > (L,L f ) (see Fig. [I]). The operator 
P acts on the full Green's functi on G (r, r') and gives the 
correlated Green's function Q T LI W^ 

g T LL , = J drdr'P(rr',TLL')G(rr'). (1) 

The integrals over r and r' are performed inside the 
sphere of size S around the correlated atom at position 
r. The subscript L can index spherical harmonics Im, 
cubic harmonics, or relativistic harmonics jm,j, depend- 
ing on the system symmetry. We always choose the basis 
which minimizes the off-diagonal elements of the corre- 
lated Green's function in order to reduce the minus-sign 
problem in Monte-Carlo impurity solvers. In general, P 
is a multidimensional tensor with one pair of indices in 
the space of local correlated orbitals (tLL 1 ) and the other 
pair in the space of the full basis set, which can be ex- 
pressed in a real space (rr') or Kohn-Sham (k, ij) basis, 
where i and j are band indices. 

The inverse process of embedding E, i.e. the mapping 
between the correlated orbitals and real-space (L,L') — > 
(r, r'), is defined by the same four-index tensor. How- 
ever, instead of integrals over real-space, its application 
is through a discrete sum over the local degrees of free- 
dom, 

E(r,r')= Yl P(r'r,rL'L)^ w (2) 

rLL'eH 

Here LL' G H means to only sum over correlated or- 
bitals. In actinides, the sum would run over 5/ orbitals, 
in lanthanides over 4/ and in transition metals over 3d 
orbitals. r runs over all atoms in the solid and r over the 
full space. Note that within the correlated Hilbert sub- 
space, the embedding and projection should give unity 
PE = I, i.e., 

J drdr'P(rr',ri 1 i 2 )P(r , r,r , i 3 L4)=^ 1 i 4 ^ 2 i 3 5rr',(3) 

while the projection from the full Hilbert space to the cor- 
related set, followed by embedding, gives the correlated 
local Green's function in real space EPG{yy') — Q(r, r') 

S(r,r')= (4) 

J2 P{t't,tL'L) j dr l dr 2 P{r l r 2 ,TLL')G{r 1 r 2 ) 

tLL'EH 



which is the central object of the functional definition of 
the DMFT described below. In general, the two opera- 
tors P and E could be different, but they must satisfy 
the condition Eq. |3|. 

The two simplest projections, namely, the projection 
on to the orbital angular momentum functions Yl, and 
the projection on to the solution of the Schrodinger equa- 
tion, can be explicitly written as 

P a {Yv',TLL') = Y L (r T )5(r-r')Y£,(r' T ) (5) 
P 1 {yt',tLL') = Y L (r T )unrr)4(r' T )Y^(r' T ) (6) 

where r r = r R T is the vector defined with the ori- 
gin placed at the atomic position i? T , and u®(r) is the 
solution of the radial Schrodinger equation for angular 
momentum I at a fixed energy E v . 

In the following, we will show that the projection P°, 
used in some implementations of DMFTE3, captures the 
full spectral weight of the correlated character L, but 
leads to non-causal DMFT equations. On the other hand 
P 1 gives causal DMFT equations, but misses some spec- 
tral weight. 

In our view, a good DFT+DMFT implementation 
should satisfy the following conditions 

(1) Correct correlated spectral weight: The projected 
density of states, computed from the projected 
Green's function, 

Pl{u) = ^[G\l{u)-GllW], (7) 

should capture the partial electronic weight in- 
side a given MT sphere at all frequencies, i.e., 

Pl{u) = pY DA (uj). In particular, Gll' must include 
the electronic weight contained in in and local or- 
bitals. Projection should not include any weigh of 
other character, nor miss correlated weight. 

(2) DMFT equations are causal: For any causal self- 
energy E, the DMFT self-consistency condition 

1 

W ~ Ei m p — S — A 

J2p^+»-h? ft -E^r 1 } (8) 

k 

should give a causal hybridization function A(w). 
Here we used projections Pk in momentum space as 
opposed to their real-space definitions in Eqs. 
and ((5), ^. 

(3) Sufficient accuracy of the hybridization function: 
The hybridization function is usually very sensi- 
tive to the choice of the projector. Therefore, we 
require that in the relevant low energy region, the 
hybridization function is similar to its DFT coun- 
terpart. Explicitly, A(w) = cj — Ei mp — (PGq)^ 1 
must be sufficiently close to its DFT estimate, 
A(oj) = uj-E imp -(P°G Q )~ 1 . HereG (r,r') stands 
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for the full Green's function G(r, r') when E = 0. 
The choice of X = is dictated by the fact that 
the hybridization A, computed by P° is not well 
behaved for S ^ 0, as we will show below. The 
motivation for using P0 in the above equation is 
that we want to project the full Hilbcrt space to 
a correlated subset with pure angular momentum, 
either / or d, but not to a mixure of characters. 

(4) Good representation of kinetic energy and elec- 
tronic density: Finally, it is crucial to faithfully 
represent the kinetic energy operator V 2 and elec- 
tronic density in real space, a feat most modern 
DFT implementations achieve. The DFT+DMFT 
implementation should not reduce the precision al- 
ready achieved in DFT underlying code. 

Downfolding to only a few low energy bands clearly vi- 
olates the condition number (3), since the hybridization 
outside the downfolded window vanishes. A more severe 
problem is that downfolding approximates the kinetic en- 
ergy operator by expressing it in a small atomic-like basis 
set, hence condition (4) is violated. Therefore, we will fo- 
cus our discussion on DFT+DMFT implemented within 
full-potential basis sets where all bands are kept at each 
stage of the calculation. Downfolding to a sufficiently 
large energy window may sometimes be helpful due to 
its conceptual simplicity, but this approach can not com- 
pute the electronic charge self-consistently, as is possible 
in our implementation. Moreover, the localized orbitals 
chosen in the downfolding procedure combined with the 
limited number of hoppings retained often cannot faith- 
fully represent the original Kohn-Sham bands. 

To be more concrete, we will give the proofs of the 
"weight loss problem" and "causality problem" within the 
full-potential LAPW basis. The equivalent derivation is 
possible for the full-potential LMTO basis. Inside the 
MT spheres, the full-potential LAPW basis functions can 
be writterPSl 



Xk +K (r) = J2AH KiL uJ K (r T )Y L (r T ) 



(9) 



where k = corresponds to the solution of the 
Schrodinger equation ui{E Vl r T ) at a fixed energy E v , 
k = 1 to the energy derivative of the same solution 
iii(E v , r T ), and n = 2, 3, . . . to a localized orbitals at ad- 
ditional linearization energies E' v , E", . . .. Here r runs 
over the atoms in the unit cell. 

The Kohn-Sham states V'ik(r) are superpositions of the 
basis functions 



^k(r) =^C i k K Xk+K(r) 



(10) 



K 



and take the following form inside the MT spheres: 

^ k (r) = YA^(k)ur(r T )Y L (i T ) (11) 



where AJ£(k) = H K 4+k l^k. or equivalently, 
/df T y L *(f T )Vik(r) - J2 K AJE(k)ur(r T ). 



The projectors ^ and ([6| can be expressed in the 
Kohn-Sham basis: 

P k (ij,rLL') = J drdr'r ik (r)P(rr',TLL')^ k (v'). (12) 

Hence, projector P° takes the form 
P£(ij,TLL') 

drdv'^(v)Y L {v T )8{r ~ r')Y£,(T' T )^ jk (r') 

= J2AIt(VAg,(k)(ur\u™) (13) 

Using projector P , we get the following expression for 
the partial density of states 

d tL ( U )= ]t ^*(k)^'(k)<uri«rV("+^-eki) 

(14) 

which exactly coincides with the DFT partial DOS. 
Hence P° satisfies the condition number (1). However, it 
does not lead to causal DMFT equations. 

To show that, consider the limit of a diverging self- 
energy, X — ¥ — ioo, as is relevant for the Mott insulators. 
Despite the diverging X, the projection must still produce 
a finite hybridization. In the case when all the bands at 
the energy of the pole are correlated, the hybridization 
should vanish. In this limit, the DMFT self-consistency 
condition (Tsl) takes the form 



(X r + A)^, 



P k (ji,rLL') 



P k (::,r'L 2 L 3 )Sl 



L3L9 



.L 2 L 3 t' 



(15) 



where :: stands for the two band indices constituting a 
matrix in ij to be inverted. Since A is finite while X 
diverges, we neglect A to obtain the condition for causal 
projection, 



Y P±(::,T'L 2 L 3 )zi 3 

L2L3T' 



(16) 



This equation must be satisfied for any matrix form of the 
self-energy S. Moreover, it has to be satisfied for each L 
and L" . We will show below that Eq. ( 16 ) is satisfied for 



a separable projection (see Eq. 19 for a definition), while 
for a non-separable projection, it likely is not. One can 
check explicitely that P° violates the condition Eq. |l6| . 
Only after applying an additional trace over LL" will the 
two matrices PS cancel. However, for any given choice of 
LL" , P° does not satisfy the causality condition. Instead 
a pole in the self-energy results in a diverging A, with 
the imaginary part having the wrong sign. The projec- 
tion P° is implemented in the qtl package^ of Wien2KP^. 
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The LDA+U implementation within Wien2K 40 also uses 
P° , but this does not cause any causality issues since the 
problem is unique to DFT+DMFT. Additionally, simple 
impurity solvers such as Hubbard-I (Ref. \T7} do not in- 
corporate a true hybridization so they also avoid issues 
with causality. 

Finally, let us mention an attractive feature of P°. 
Within this scheme, the self-energy is independent of the 
radial distance from the atom r T , having only angular 
dependence in the form S(r,r'). This matches the con- 
ceptual fact that the impurity solver within the DMFT 
framework can not determine the radial dependence of 
the self-energy. The impurity solver can only be used to 
obtain the angular dependence of £ by determining the 
expansion coefficeints Sll'- ln the absence of any knowl- 
edge of the radial dependence of E, the natural choice 
is a constant function, independent of radius r r . Since 
£(r,r') is a function of two vectors, a radial delta func- 
tion would be an obvious choice. However, issues with 
causality preclude the use of this projection. 

The second projection P 1 of Eq. |6| takes the following 
form in the Kohn-Sham basis: 

P£(ij,rLL') = 

^AT£(k)Ag,*(k)(ur\uf)(4\uJ, K '). (17) 

The partial density of states computed from the corre- 
lated Green's function using P 1 is 

D tL {lo) = (18) 
£ ^(k)A[«'*(k)«"| U 0)(^K»V( W + M - e kl ) 



Comparing Eq. (18) with (14 1, we notice that (uJ K \u TK 



is replaced by 
spectral weight 



i \ u i ) 

which leads to incorrect 
the original 



In particular, for k = 1 
overlap in Eq. (14 1 is (uj\uj), while the overlap obtained 
by P 1 , vanishes. 



Causality is not violated for any projection P, which 
is separable, i.e., can be cast into the form 



P\i Jl TLL') = U%Ufl; 



(19) 



The condition Eq. (16) can then be expressed as 



1 = ?7 krt (C/ kr S r t/ kTt )" 1 C/ kr S T 

k 



(20) 



which is clearly satisfied when C/ kr is invertible ma- 
trix because W{UEW)- X UE WUiWU)' 1 = 1. This is 
satisfied when the Kohn-Sham Hilbert space is of larger 
dimension than the correlated Hilbert space. The projec- 
tion P 1 leads to causal DMFT equations, and therefore is 
a better choice than P°. However, some spectral weight is 
lost at energies away from the linearization energy E v . To 
this end, we also implemented an alternative projection 
within Wien2K package 31 , which preserves both causal- 
ity and spectral weight. This projector is given by 



P\vr\rLL') = ]T ^k(r)^(k)( U r|< )(<°|</)^;*(k)^ k (r') x 



E 



K\K2 iL 



\ \ V A TKl A TK 



TK 2 * /,,TKl |„,T0 



r°><« 




ATKl* ATKl . 
ATK±* A TK -2I 

k 2 •A-jL' jL' < 



5 >K 



(21) 



Here index L runs over the local basis in which the green's 
function is minimally off-diagonal (cubic harmonics or 
relativistic harmonics) . 



U. 



with 



iL \l ATKl ATK 2 * l-.TKi I „,t0\ /,,tOL,™. 2 \ 



The projector is separable, as postulated in Eq. (19), 
and the transformation U is 



(22) 



(23) 



E K1 , 2 ^f^r<< K >r°><<K K2 > 

Hence the DMFT equations are causal. Moreover, 



P£(ii,LL) is identical to P^(ii,LL) and hence the par- 
tial density of states Dl(uj), obtained by P 2 , is identical 
to Eq. ( |14[ ). Hence the projection correctly captures the 
partial spectral weight. Knowledgeable reader would no- 
tice that the projection is slightly non-local because SJ L 
is weakly momentum dependent. At energies where u 
or local orbital substantially contribute to the spectral 
weight (away from the Fermi level), we give up locality 
in expense of correctly capturing the spectral weight. 

All projection schemes lead to slightly non- 
orthonormal correlated Green's function. This is 
because the interstitial weight is not taken into account 
and because the full potential basis is overcomplete. To 



G 



have an orthonormal impurity problem, we compute 
the overlap '}2iiP 2 {'ihTLL') = O t ll , and renormalize 
P 2 {i^tLL') -> E ili2 (^)M 1 P 2 (ij,rL 1 L 2 )(^) i2 ^. 

Finally, we remark that the segment of our code which 
builds projections P°, P 1 and P 2 within Wicn2K 31 is 
based on the qtl package of Pavel NovalP^. 

Similar projections within LDA+DMFT method were 
proposed before. In particular the method by B. Amadou 
et.al.^ proposed to construct the Wannier functions for 
the correlated subset only, while the DMFT equations 
were solved in the Kohn-Sham basis, restricted to some 
subset of low energy bands. The local orbitals used 
for the projection were either all-electron atomic partial 
waves in the PAW framework, or pseudo-atomic wave 
functions in mixed-basis pseudopotential code. Hence, 
in the language of projectors, the method was similar to 
choosing the projector to be P = \xSm ><: XkmL where 
Xkm 1S the the partial waves or pseudo-atomic wave func- 
tion. While this method is clearly causal, it looses spec- 
tral weight of the correlated angular momentum charac- 
ter. Moreover, the implementation of the method did 
not allow the self-consistent evaluation of the electronic 
charge. The method of Anisimov ei.aZPSl also proposed a 



construction of the Wannier functions using an arbitrary 
set of localized orbitals. In their work, the LDA Hamil- 
tonian was truncated to Wannier representation for the 
purpose of obtaining the DMFT self-energy. This sim- 
plifies the self-consistent DMFT problem, but makes it 
impossible to implement the charge self-consistency. Fi- 
nally, Savrasov et.alP^proposed a projector particular to 
LMTO basis set, for which causality was not proven. 



III. DFT+DMFT FORMALISM 

To derive the DFT+DMFT equations, we define a 
functional of the correlated Green's function Q(r,r') and 
extremise it. The correlated Green's function Q(r,r') is 
defined by Eq. Q, and the functional to be extremise is 

T[G,p] = -Trln(G _1 ) - Tr[£ tot G] + $[G,p\, (24) 

where Tr runs over all space (orbitals, momenta) and time 
(frequency) . The quantities apprearing in the above func- 
tional are 



Gz\r, r') = [uj + M + V 2 - V^r)] S(r - r') - (r, r') (25) 

E£"(r, r') = [V H (r) + V xc (r)} 6(r - r') + [E„(r, r') - E DC 6(r - r')] 0(r < S) (26) 

$[G,p] = $u[p] + ®xc[p] + <S>dmft[9] - ^dc[G] (27) 
P = Tr[G] 



where Tr is trace over time only (not space), V ext is the 
potentials due to ions, Vh,Vxc ar e the Hartree, and 
exchange-correlation potential, respectively. $dmft[G] 
is the sum of all local two particle irreducible skeleton di- 
agrams constructed from Q, and the Coulomb repulsion 



U (screened by orbitals not contained in Q), and &dc is 
the double counting functional. 

We assume that the Coulomb interaction U has the 
same form as in the atom, i.e., 



U 



^ ^ 21: 

L a ,..L t i,,m,er<7 / k—0 



^ 2 k + \ ( Y LA Y krn\YLj(Y L jY£ m \Y Ld )fl aa fl bCr ,f Lda ,f LcC 



(28) 



however, the Slater integrals are reduced due to screening 
effects. Typically, we renormalize F 2 ■ ■ ■ F 6 by 30%, from 
their atomic values, while F°, being renormalized more, 
can be estimate by constraint LDA or constraint RPAHD. 



To extremize the functional Eq. (24), we take Q and p 



as independent variables, and use the following functional 
dependence: E[£], $dmft[G], E dc [G], $dc[G] are func- 
tionals of Q. Consequently, G is also a functional of Q, 



i.e., G[E[£]]. On the other hand, V H [p], V xc [p], <f> H \p], 
&xc[p] are functionals of the total electron density, hence 
G is also a functional of p since G\Vh[p] + V xc [p]]. Finally 
it is easy to check that 



Tr[E tot G] = Tr[(V H + V xc )p] + Tr[(E - E DC )G]. 
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With the above functional dependence in mind, mini- 
mization with respect to Q gives 



S — Edc — 



S$dmft[G] 6$ dc [G] 



sg 



sg 



and minimization with respect to p leads to 



V H + V xc = 



dp 



dp 



Hence the Hartree and exchange-correlation potential are 
computed in the same way as in DFT method (note how- 
ever p is electron density in the presence of DMFT self- 
energy) , while the DMFT self-energy is the sum of all lo- 
cal Fcynman diagrams, constructed from g and Coulomb 
interaction U. 



To sum up all local diagrams, constructed from g and 
screened Coulomb interaction U, we solve an auxiliary 
quantum impurity problem, which has g — Gi mp as the 
impurity green's function, and E as the impurity self- 



energy E 
G 



imp 



imp 



= l/(» 



E 

- E, 



The impurity Green's function is 
A), hence the DMFT 



■imp 



self-consistency condition reads 



P(iuj + p-H 



DFT 



-.BE))- 1 = (iu-Et, 



A)- 



(29) 

where E = E — Edc> an d -Edc is the interaction in- 
cluded in DFT (double counting). The self-consistency 
condition takes the explicit form 



f drdr'P(rr', tLL') J [iu + p + V 2 - V KS (r)] S(r - r') - V P(r'r, tL x L 2 )^ Li 1 

J (ry)<Sr { LlL2eH J 

= [{iw-EJ^-^-A^- 1 



r 



(30) 



where Vks = V ex t + Vh + V xc and S is the muffin-tin 
radius. 

For efficient evaluation of the DMFT self-consistency 
condition Eq. (30), we choose to work in the Kohn-Sham 



(KS) basis. At each DFT+DMFT iteration, we first solve 
the KS-eigenvalue problem 



[-V 2 + V KS (r)] Vki(r) = e k iVki- 



(31) 



Then we express the projection P in KS basis, 
Pk(ij,rLL'), where i,j run over all bands. We then per- 
form the embedding of the self-energy, i.e., transforming 
it from DMFT base to the KS base 



Sk,ij(w)= ^2 P kT (ji,TL 2 L 1 ) Z T Lil 



US) (32) 



In KS-base, we can invert the Green's function Eq. (30), 



to obtain the practical form of the self-consistency con- 
dition 



kij 



LL' 



1 



iui — EJ — E T (o;) — A T (w 



iu + p — e k — E k (w)) 1 [33) 

Jji 

(34) 



LL' 



This is of course equivalent to Eq. (30). Finally we solve 



this self-consistency equation for a given self-energy E(w) 
to obtain the hybridization function A T and the impurity 
levels EJ mp . 



We note in passing that the self-energy E(w) is a 
complex function, and its imaginary part is related to 
the electron-electron scattering rate, which is very large 
in correlated materials. In Mott insulators, it is even 
diverging. Hence the DMFT "effective Hamiltonian" 
£k + Ek(w) can not be diagonalized by standard methods 
to obatin a set of eigenvalues, i.e., bands. The eigen- 
values are complex and hence only the spectral weight 
A(k,w) = (G k H - G k (w)/(27ri) is a well defined quan- 
tity. The absence of well defined bands in correlated ma- 
terials makes computational techniques more challeng- 
ing. For example, the calculation of the chemical poten- 
tial is far more demanding because one can not assign 
a unity of charge to each fully occupied band. Rather 
all complex eigenvalues, even those which are far from 
the Fermi level, need to be carefully considered. This 
point will be addressed below in section [TV] item 5. Fur- 
ther, the tetrahedron method^, a very useful technique 
to reduce the number of necessary momentum points in 
practical calculation, is not applicable since it needs real 
eigenvalues. We address the necessary generalization of 
this method is chapter [V] 

Note that generalization of the projector and the 
LDA+DMFT formalism to cluster-DMFT is very 
straightforward. One needs to increase the unit cell to in- 
clude more sites of the same atom type. The self-energy 
and the Green's function become matrices in index r, 
i.e., E£^,, QWi- The transformation P is also straight- 
forwardly generalized to matrix form P^ij; tLt'L'). The 
only difference in the definition of the projector Eq. (21 ) 



is that A T L < is replaced by A T L , (A T L remains unchanged), 
which amounts to the integral over two different spheres 
around two atoms of the same type. Finally, in cluster- 
DMFT case, the self-energy in KS-basis Eq. Q has 
to be summed over both t and r', and self-consistency 
condition Eq. (34) becomes a matrix equation in t,t'. 



The challenging part of the cluster-DMFT formalism is 
in solving the cluster-impurity problem. In combination 
with impurity solvers based on the hybridization expan- 
sion (discussed below) the computational effort grows ex- 
ponentially with the number of correlated sites. In the 
weak coupling impurity solvers, the computational effort 
grows as a power-law, however, these techniques usually 
can not reach the interesting regime of strong correlations 
and low temperatures. 

The major bottleneck in evaluating the DMFT self- 
consistency condition in our method is the multiplication 



of the projector Pk T (ij, LL') with E in Eq. (32 1 and mul- 



tiplication of projection with Green's function Gkjj in 
Eq. (33). Since projection P 2 is separable, one can write 
the operation in terms of matrix products. Still, these 
sums run over all k-points (typically few thousands) and 
all frequency points (typically few hundreds). 
For the efficient implementation of the set of Eqs 



(32) 



and ( 33 1 , we first notice that the transformation P (or its 



separable part U) is very large and is not desirable to be 
written to the computer hard disc. Hence we generate it 
only for one k-point at a time, and evaluate both prod- 
ucts at this particular k-point. Non- negligible amount of 
time is necessary to generate the transformation Eq. (21 ), 



and because this transformation does not depend on fre- 
quency, it needs to be used for all frequencies in Eqs. ( 32 ) 
and ( 33 1 . Hence paralization over frequency is not im- 



plemented, while paralization over k-points is. 

Note that because of the sum over atoms (t) in 
Eq. (32), the transformation for all atoms needs to be 
computed first, and only then the sum in Eq. ( |32[ ) can 
be evaluated and the self-consistency condition Eq. (34) 
can be inverted. 

To optimize the sum in Eqs. ( 32 ) and ( 33 ), one can no- 



tice that local quantities like self-energy and local green's 
function possess a large degree of symmetry when written 
in proper basis (real harmonics, relativistic harmonics): 
many off-diagonal matrix elements vanish, and many ma- 
trix elements are equivalent. For example, in a d sys- 
tem with cubic symmetry, one has only two types of self- 
energy t 2g and e g . Hence, instead of summing over 10 x 10 
matrix elements in Eq. ( |32[ ) , one can rewrite the sum over 
two matrix elements t = (0, 1), i.e., 



(35) 



where P kT (ji,t) = Es(L 1: L 2 )=E(t) PkAji,L 2 L 1 ) and the 
indices Li,L 2 here stand for the real harmonics rather 
than spheric harmonics. The later transformation is in- 



dependent of frequency, while the sum Eq. ( 35 ) needs to 



be performed for all frequencies, hence the compact form 
of the transformation saves a lot of computer time. 



IV. THE ALGORITHM 

The implementation of the DFT+DMFT algorithm is 
done in the following few steps: 

1) p(r): We converge the LDA/GGA equations to get 
the starting electronic charge p(r). We use the non- 
spin polarized solution as starting point. In the 
ordered state, the DMFT self-energy is allowed to 
break the symmetry, while typically the exchange- 
correlation potential is not allowed to break the 
symmetry (LDA rather than LSD A). 

In this preparation step we also obtain good esti- 
mates for the Coulomb repulsion U (which is rep- 
resented by Slater integrals F°, F 2 , F 4 and F 6 ). 
Slater integrals are computed by the atomic physics 
program of Ref. 2H and they are scaled down by 
30% to account for the screening in the solid. The 
F° terms is very different from the atomic F° and 
is obtained by constraint LDA calculation, or con- 
straint RPA calculation^. 

2) ^ki(r): We solve the DFT KS-eigenvalue problem 

(-V 2 + V KS (r))i; ki (r) = Mr)eg FT 

to obtaine KS eigenvectors, core, and semicore 
charge, and linearization energies E v . 

3) Ell' : We start with a guess for the lattice self- 
energy correction E(w) = E(w) + Eqo — Ed c (here 
E is the dynamic part of the self-energy with the 
property E(oo) =0). A reasonable starting point 
is E(w) = and Ed c = (£>oo)- The potential in the 
first DMFT iteration is thus the DFT potential. 

4) E k ,ij: Next we embed the DMFT self-energy 
T, ll ,(oj) (shifted by double counting) to Kohn- 



Sham base by the transformation Eq. ( 32 ) to obtain 

^k,ij(w). 

5) /x: Using the current DMFT self-energy E(w), and 
the current DFT KS-potential VkSi we compute 
the current chemical potential. This is done in the 
followin steps: 

— Complex eigenvalues £kz(w) of the full Green's 
function are found in the large enough energy 
interval (at least [— 2U,2U]) by solving 



5 kii (w)]C)5(a;)=C?f( t( ;) ekI (w). 



Here Cji are DMFT eigenvectors expressed in 
KS base. The DMFT eigenvalues outside this 
interval are set to DFT eigenvalues. We need 
only eigenvalues in this step, but not eigenvec- 
tors. 
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The chemical potential is determined using 
precomputed complex and frequency depen- 
dent eigenvalues Eki,^- On imaginary axis we 
solve 



N„ 



i 



and on real axis we solve 
I ui 4- 



iV, 



w + m - £k;(w) 



If enough k-points can be afforded, we use 
special point method, otherwise the "complex 
tetrahedron method" can be used (see chap- 
ter 0. 

For numerical evaluation of the real axis den- 
sity, we discretize the integral 



N„ 



1 



im^r/ME 



k( 



with ai = (uji + Wj_i)/2 and bi = 
(u>i+i + u>i)/2. When using the special point 
method, the integral over frequency is evalu- 
ated analytically, and the terms of the form 
log(di + \i — £kz(wi)) are summed up. Alter- 
atively, we sometimes use the complex tetra- 
hedron method, where the four-dimensional 
integral is evaluated analytically (see chap- 
ter |v} 

When DMFT is done on imaginary axis (using 
imaginary time impurity solvers), we evaluate 



N 



f(4i - m) + 



ki 



2T x: 

0<ui n <uj N kZ 

1 



y\ t t ] 



/ £ ti - A* 
- arctan — 

7T V LO N 



1 



ai'rtai) I — H) 

7T \ uj n ■ 



Here e^/ is the real part of the eigenvalue at 
arbitrary frequency. We choose the lowest or 
the last Mastubara point. Again, the tetrahe- 
dron method can be used for momentum sum. 

For Mott insulators, the above described method is 
not very efficient, because even a small numerical 
error in computing N va i places chemical potential 
at the edge of the Hubbard band, either upper or 
lower. This instability usually does not allow one to 
reach a stable self-consistent solution. We devised 
the following method to remove this instability: 

— The diagonal components of the self-energy 
were fitted by a pole-like expression S^l' = 



5joo + 



W L 



— Next, we neglected broadening of the pole 
(r^), which should be small in the Mott in- 
sulating state. We computed a quasiparticle 
approximation for the Green's function G k p , 



DFT_y .._r/ kT 



1 



is part 



Pi v ■ ■ ±j — j L 
L J 

(37) 

of the projector 



where 
P^(ij,tLL' 

— The above Green's function formulae can be 
cast into a block form 



UflU^l* defined above. 



Gl P = 



= {iu,-H?). 



j^ + Eoo £/ k VVF 
v / VFf7 krt P 

Here H^ p is the quasiparticle Hamiltonian 
which can be diagonalized to obtain the quasi- 
particle bands. We notice that the num- 
ber of quasiparticle bands of the Mott insu- 
lator is larger then the number of Kohn-Sham 
bands because Mott insulators have at least 
two Hubbard bands. The quasiparticle bands 
are not very accurate away from the Fermi 
level, however they are sufficiently acurate at 
low energy and allow one to identify gaps at 
the Fermi level. Once a gap in the spectra 
of is identified, the charge is computed 
using the full DMFT density matrix to verify 
the neutrality of the solid. If the solid is neu- 
tral when chemical potential is in the gap, the 
chemical potential is set to the middle of the 
gap. 

6) A: Impurity hybridization function A(w) and im- 
purity levels Ei mp are computed in this step. 



We use equation (|33| to get Q L lt and we use the 



high frequency expansion of both equations ( 33 ) 
and (34) to determin impurity levels 



E 



imp LL , 



= -E 



DC<JLL' 



DFT 
ki 



kf 



7) Ei mp : Impurity solver uses Ax,zy(u;), E, 



■imp j 



and 

Coulomb repulsion U (which is represented by 
Slater integrals F°, F 2 , F 4 and F 6 ) as the input 
and gives the new self-energy Eii'(w) as the out- 
put. 

Currently we integrated the following impurity 
solvers: OCA (see chapter VII B), Non-crossing 
approximation (NCA), Continuous time quantum 
Monte Carlo (CTQMC^. The latter is imple- 
mented on imaginary axis, and the former two on 
real axis. 

Before the impurity solver is run, we exactly diag- 
onalize the atomic problem in the presence of crys- 
tal fields, to obtain all atomic energies E m and the 
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matrix elements of electron creation operator in the 
atomic basis (m\f' a \n). Since the impurity levels 
can change during the iteration, the crystal field 
of the atomic problem can change as well. In case 
of /-systems, the crystal field splittings are small 
and one can assume that they do not change sub- 
stantially from their DFT value. Hence the exact 
diagonalization can be done only once at the be- 
ginning. For the ci-systems, the crystal field split- 
tings are larger, and this approximation is in gen- 
eral not necessary satisfactory, hence the exact di- 
agonalization needs to be repeated in the charge 
self-consistent cycle. A special care needs to be 
taken here when using CTQMC. To speed up the 
convergence of CTQMC solver, we typically start 
simulation with the status of the kink distribution 
from previous DMFT step. Since exact diagonal- 
ization can reorder eigenstates, these kinks need to 
be properly renumbered, to efficiently restart sim- 
ulation. 

Soo: It is very hard to achieve reasonably pre- 
cise self-energy at high frequency with impurity 
solvers based on hybridization expansion. How- 
ever, to correctly compute electronic charge, it is 
crucial that the self-energy at high frequency ap- 
proaches its Hartree-Fock value and the impurity 
Green's function and self-energy at large frequency 
properly behave. Hence we correct Eoo at each it- 
eration. This is quite straighforward, given the fact 
that impurity solvers determine the impurity den- 
sity very precisely. This steps only corrects the high 
energy tails of the impurity green's function and 
impurity self-energy, while we make sure that the 
low energy part, which is computed very precisely 
by these methods, is not altered. 

In the case of CTQMC solver, we compute the 
atomic Green's function using CTQMC probabil- 
ities for each atomic state (see Ref . [2H for details) . 
The high-frequency tails of the self-energy can then 
be computed. These analytic tails are then used in- 
stead of noisy QMC data. 

In OCA and NCA impurity solvers, we project out 
very high excited atomic states. This has negligible 
effect on the low energy physics, however, it results 
in a missing weight at high frequency, and hence 
wrong self-energy at infinity. To correct for this 
deficiency, we add two lorentzians to the impurity 
Green's function 



0M = 



A(x)dx a\ a,2 

LU — X UJ — 6\ + iT LJ — 62 + iT 



typically with e% < —U and t2 > U. Here we omit- 
ted the subscript LL' for the impurity Green's func- 
tion Gll> for clarity. The parameters ai, 02,61,62 
are determined by the following constraints: 



normalization: m + a% 
the integral of A(x). 



a 2 = 1 , where m is 



— density: n + a\ — n exactl where n 

I A(x)f(x)dx and n exac t is the impurity den- 
sity determined by the impurity solver in an 
alternative, more precise way (from pseudo- 
particle density). 

- Soo: mi + aiei + a 2 e 2 = E imp + Soo, where 
mi is the first moment mi = J xA(x)dx. 

Once the following three constrains are satisfied, 
the self-energy at high frequency approaches its 
Hartree-Fock value, and the spectral function re- 
spects the total impurity density. 

9) S: Using the new impurity self-energy, we deter- 
mine the new lattice self-energy S(w) = £(w) + 
^-Edc, where E DC = U(n-\/2)- J(n/2-l/2), 
with n the correlated nominal occupancy. 

10) goto 4 '■ If the convergence of charge is hard to 
achieve, we iterate the DMFT loop a few times. 
We call this loop the DMFT loop. If the DMFT 
loop is to be iterated, jump to 4. 

11) [i, p(r): The eigevalue problem is solved for all mo- 
mentum and frequency points, 



3 



.DFT, 



Here we evaluate both, eigenvalues and eigenvec- 
tors. Since this is a non-hermitian eigenvalue prob- 
lem, the left and right eigenvectors are not complex 
conjugates of each other. We use notation C£ff for 
the right and CQh for the left eigenvector. 

Using the DMFT eigenvalues, we recompute the 
chemical potential as in 5. 

We then recompute the electronic charge from the 
DMFT eigenvectors 



^kiw(r) 



E 



where ip-^i are Kohn-Sham eigenvectors (solutions 
of the LDA eigenvalue problem). The electronic 
valence charge on real axis is 



p va i(r) = -hm^2 I V>kL(r) 

kZ J 



fiUj)dUJ -^(r) 



w + n - 



and on imaginary axis is 
p val (v) =TJ2 ^kk(r) 1 



We compute the electronic charge using similar 
technique as used above to compute the chemical 
potential. The electronic charge is 



ft,«rf(r)=2^ ki (r)^ J .(r)Wgj 



DMFT 
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The weights W-^^ FT on real axis are combuted as 



and evaluated by 



DM FT _ p <J p ]i r> u p i . 



k.,., — / , Ckil ( 'kjl "'kip 

Ip 



with 



1 f p 1 
Wki p = f(uj p )lm / du — ■ 

and a p = (w p + w p _i)/2, 6 p = (w p +i + w p )/2. 

On imaginary axis we evaluate the weights by the 
following expression 



DMFT 



. v ?cj„ + /it - e k iw„ w*>n + M _£ Ww / 

+ E C k«*/( £k ^ -M) 

i 

Note that the DMFT density matrix W£j$ FT is 
a hermitian matrix in Kohn-Sham band indeces i 
and j. Hence, we can use eigenvalue techniques for 
hermitian matrices to decompose W into 



i 



The LDA+DMFT electronic charge can then be 
evaluated by rotated Kohn-Sham vectors, and 
DMFT weights w k j by 



j 



Hence, the code to compute the LDA charge can 
be simply converted to compute the DMFT charge 
by just replacing the Kohn-Sham LDA weight by 
DMFT weight and by rotating the Kohn- 

Sham eigenvectors by the above computed eigen- 
vectors Uk. 

Finally, the DFT core and DFT semicore charge 
is added to the valence charge, and the resulting 
total charge is renormalized in the standard way, 
such that the charge neutrality is satisfied to high 
accuracy. 

12) i? tot :The total energy is computed on the output 
density p(r), using the low temperature limit of the 
functional Eq. Q evaluated on the DFT+DMFT 
solution: 



E, 



total 



Tr[(-V 2 + V ext )G] + ^Tr [EG] +E H + E XD -& DC 

For computation, the formula is cast into the fol- 
lowing form 

Etotai = Tr[(-V 2 + V KS )G] - J (V ff (r) + V xc {v))p{v)dv 
+ E H + E xc + ^Tr[SG] - $ DC 



rr — c DFT w 

Etot - 2^ e k' W \ 



DMFT 
k,ii 



(V H (r) + V xc (r))p(r)dr 



+ E H + E xc + E l ™f enUal - <5> DC 

where Wg MFT is the DMFT density matrix de- 
fined above, and 



potential 2 



\t ^lIMqPm (39) 



,tLL' 



is the impurity potential energy which can be com- 
puted very precisely by most impurity solvers, such 
as CTQMC or OCA. For example, in CTQMC we 
sample probability for each atomis state P m . Us- 

J potential 



ing these probabilities, we can evaluate E 1 ™ 11 " 



Ep rpatom 



LL> 



imp ^jmp 



LL 1 ^LL' L'L 



13) mix: The total electronic charge is mixed with the 
charge from previous iterations using multi-secant 
mixing of Marks and LukeP. 

14) DFT: In this step, we recompute the DFT poten- 
tial (hartree, exchange-correlation potential), the 
Kohn-Sham orbitals and linearization energies. 

15) goto 11: If the self consistency is hard to achieve, 
jump to 11 and determine the best electronic charge 
p(r) on the current impurity self-energy S. We call 
this loop the LDA loop. 

16) goto 6 If the electronic charge and self-energy are 
not converged, jump to 6. We call this loop the 
charge loop. 



V. COMPLEX TETRAHEDRON METHOD 

The calculation of the electronic density, as well as 
the correlated Green's function, requires precise eval- 
uation of integrals, which contain diverging poles. In 
systems with many atoms per unit cell, one can not af- 
ford enough k-points to get hybridization function A(w) 
smooth on a scale of temperature T without introduc- 
ing artifical broadening larger than T. Hence, to avoid 
artifical broadening larger than the low energy scale, we 
need to use alternative summation over momentum. The 
tetrahedron methocP^ is used in this case. In the con- 
text of DFT+DMFT, an aditional complication is that 
the eigenvalues are complex numbers. Although the ana- 
lytic formulas for the integration over a tetrahedron can 
straighforwardly be evaluated, and are given in appendix 
A, a more severe problem is the interpolation of the mul- 
tidimensional complex functions in momentum space. 
Below we give details on a method to overcome this dif- 
ficulty. 
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Computation of the Green's function requires the eval- 
uation of the following integral 



Ca 



which can be rewriten as 



where the sum runs over all tetrahedrons t, and integral 
needs to be performed over the particular tetrahedron 
t. i is the band index. The linear interpolation of Cjk 
and linear interpolation of in momentum space leads 
to analytic formulas for the weight functions w(i,\s.,cj) 
(given in appendix A), which can be used to evaluate g 
to higher precision by g = ^ k w(i, k, cj)Cik- 
Similarly, the electron density is computed by 



N„ 



dojfjuj) 
u> + fi- e iku ' 



We take a frequency mesh, which is sufficiently dense at 
zero freqeuncy that it can resolve the fermi function f(co), 
and we approximate 



N val = --ImV/(^) f d 3 k f 



duj 



It J(u j +u j - 1 )/2 W + M^ e ikw 

ilm^/K0^(k,^±^,^^ 



140) 



Here the integral J. is the integral over a particular tetra- 
hedron i. The weights can again be computed analyti- 
cally and are give in Appendix A. 

To evaluate the integral over a tetrahedron t, which 
has corners in momentum points ^1,^2,^3,^4, we need 
to interpolate the eigenvalues e ilkl , £i 2 k 2 , £i 3 k 3 , ei 4 fc 4 inside 
the volume of the tetrahedron. Since there are many 
crossing bands (index i), it is not at all simple to find a 
good interpolation of e ik inside the tetrahedron. 

In the standard tetrahedron method, where eigen- 
values are real numbers, one sorts the eigenvalues at 
each fc-point, to get the vector of increasing energies 
£ i,fc, £2,fc, • • • , and then one linearly interpolates each 
sorted component of the vector ei fkl , ei,fe 2 , ej,fe 3 , £i,/c 4 in- 
side the tetrahedron. Hence all crossings are avoided. It 
is however important that no artifical crossings are ob- 
tained in the interpolation, because a crossing gives a 
diverging contribution to the integral. 

Complex eigenvalues, which appear in DFT+DMFT, 
can not be sorted. Hence the interpolation is not at all 
simple. A reasonable attempt would be to sort eigen- 
values according to their real parts, and just neglect 
their imaginary parts when sorting. It turns out that 
in strongly correlated regime, where the self-energy be- 
comes very large at some frequency points, the error in 
tetrahedron method is so large that the hybridization 



function can become non-causal in such points. Due to 
this non-adequate interpolation, the Green's function has 
a lot of noise, superimposed on a smooth curve. How- 
ever, hybridization function, which is many times more 
sensitive than the Green's function, has unbearable large 
error, which cause enormous error in the solution of the 
impurity problem. 

To overcome this problem, we implemented a special 
type of smooth interpolation, based on the idea that the 
absolute value of the energy should not change much from 
one k-point to its neighboring k-point. For each tetrahe- 
dron, we minimize the following functional 



i (ot,0)Epairs 



— mm 



(41) 



where the 6 pairs of the tetrahedron corners are: 
(1, 2), (1, 3) • • • (3, 4), and i runs over all bands. We min- 
imize the functional with respect to the order of eigen- 
values in all corners of the tetrahedra. 

To minimize the above functional, we can choose an 
arbitrary order of bands in the first fc-point k\ , and then 
we have to permute the components of the other three k- 
points {k2,kz,ki). Hence the number of all possible trial 
steps is (n!) 3 , where n is the number of bands, and is 
typically of the order of few hundred. Obviously, not 
all arrangements of the eigenvalues can be tried. Our 
algorithm for sorting the eigenvalues is 

1 Sort the eigenvalues according to their real parts. 

2 Use Metropolis Monte Carlo method (for T = 0) 
to flip components of a vectors e k ,i e k ,j. Try 
to flip components in any of the momentum points 
k 2 , ■ ■ ■ , fc 4 . 

The trial steps are chosen in such a way that the probabil- 
ity for flipping two eigenvalues, which have very different 
real parts, is very small. We typically choose an expo- 
nential distribution function with probability P(\i—j\) oc 
exp(|i- j|/5). 



VI. TRANSPORT CALCULATION USING 
DFT+DMFT 

In this section, we will give the efficient algorithm 
to compute the DC conductivity within DFT+DMFT. 
The higher order transport coefficients can be computed 
along the similar lines, although the computation be- 
comes more technically involved. 

The DC-conductivity can in general be expressed by 



\1V 



^^^(w + iS) 



(42) 



where the current-current correlation function \ IS ex_ 
pressed diagrammatically through the electron Green's 
functions and the current vertex function by 
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-T 



E 



v 



k/j s~tPiPi i 
P1P2 kir 



GZ («OGk - ^njr^ ( k ^ w «) 



(43) 



kav m ,pi,p 2 ,p' 1 ,p l 2 

Here r(kz4 n ,u; n ) is the current vertex function, which satisfies the integral equation 

■ k " - ^ ;(b mi kV>„)^(i 1 /J^(i i / m - io;„)r^ 4 (k'^,a;„)(44) 



jcra 

P'lP^P'sP'i^ 



r 



and J(ki/ m , k'z/^; w n ) is the particle hole irreducible ver- 
tex, whose limit at zero frequency and Fermi momenta is 



the Landau interaction function, 
by 



iv are velocities, given 



P1P2 



2m 



(0k Pl |V 1/ |V>kp 2 



All quantities are expressed in a Bloch-basis, for example 
the Kohn-Sham basis, which diagonalizes the static part 
of the action. 

In general, the two particle vertex function is very dif- 
ficult to compute. In some cases, the vertex corrections 
vanish and the transport quantities can be computed 
from the lowest order bubble diagram. 

If self-energy is momentum independent, and the single 
band approximation is appropriate, the vertex correction 
vanish, as shown by Khurana^. In multiband system, 
the following set of conditions arc sufficient for the vertex 
correction to vanish: 

1) The irreducible vertex function is local, i.e., 
/(kf, kV;u; n ) does not depend on k or k'. 

2) Velocities are odd functions of momentum, i.e., 

3) Green's function is even functions of momentum, 
i.e., G_ k = G k . 



Under the above conditions, it is clear from Eq. (44 1 that 



only the zeroth order term remains and vertex is un- 



Dynamical Mean Field Theory the self-energy operator is 
approximated by a purely local quantity. However, the 
local approximation is made in a localized basis. The 
self-energy in the Kohn-Sham basis is given by Eq. (32 1, 
and is obviously momentum dependent. In general case, 
the resulting self-energy £k+-PkrE is not an even function 
of momentum, and hence G_k 7^ Gk. 

Due to difficulties in computing the two particle vertex 
function to high accuracy on real axis, the vast majority 
of theoretical calculations ignore the vertex corrections 
to conductivity. At present it is not clear how impor- 
tant the vertex corrections to optical conductivity and 
transport are in correlated electron materials. They are 
likely small because they vanish at low energy, where 
an effective single band approximation is possible. And 
they are also small at intermediate energies where the 
interband transitions give major contribution to optical 
conductivity. However, a thorough investigation of the 
vertex corrections and consequently appearance of exci- 
tons in correlated materials is a very interesting avenue 
for future research. 

In the absence of vertex corrections, the current- 
current corelation function Eq. ( 43 1 becomes 



v ° k J 

Tr (p k (y)v^p k (y-Lj)v^) (45) 

where pk = (G^ — G)/(2iri). Both spectral density pu 
and velocity are matrices in orbital indices and trace 



k . Consider the first order term is taken over the orbitals and spins in Eq. (45). Finally, 



renormalized T(k) 

Sk' I Gk'Gk'« k in Eq. (|44|) or the second order term 



the real part of the DC conductivity is given by 



Ek'k» J GvGvI G k »G k »v k " in Eq. ( 44 1 . The func- 



tion being summed is odd in k' and k", respectively, and 
hence the terms vanish. 

Under which circumstances the above three conditions 
are met? The first condition is exact in the limit of infi- 
nite dimensions. Thus in Dynamical Mean Field Theory, 
the irreducible vertex is local. For many three dimen- 
sional systems, it is believed to be an excellent approx- 
imation. However, the velocities are not necessary odd 
functions of momentum, in particular, they are obviously 
nonzero in strict atomic limit, thus violating the condi- 
tion (2). Finally, the third condition is obviously satis- 
fied in single band theories with local self-energy, where 
Gk(uj) = l/(w + p — ek — E(w)) because e_k = £k- In 



V ° k 



dy 



dy 



Tr( p k (y)v^p k (y)v^)m 



The dynamic self-energy is computed by an impurity 
solver, which is implemented either on the real or imag- 
inary axis. The most precise impurity solvers, such as 
CTQMC, are implemented on imaginary axis, hence we 
would like to formulate the method also for the case of 
imaginary axis self-energy. Since the DC transport is 
sensitive to the behaviour of the self-energy at low fre- 
quency, we take the power expansion for E(iw) and we 
determine the coefficients directly on imaginary axis 



E(w) = E(0) + (1 - Z- X )u - iu?B + 



(47) 
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For the DC conductivity, the expansion to the quadratic 
order is quite accurate. However, for the thermoelectric 
power, the truncation at quadratic order is not sufficient 
since the qubic terms in the self-energy expansion (the 
asymmetry of the scattering rate) is crucial even at low 
temperature (see RefP^). 

We first embed the quasiparticle renormalization am- 
plitude Z and scattering rate B to the Kohn-Sham basis 
using Eq. Q, i.e., Z^ 1 = P kT Z- x and B k = P kT B. 
Then we can express the low energy electron Green's 
function in the Kohn-Sham basis as 



G k (w) = (wZ k 1 + (j, - £(0) - e k + iuj 2 B k y 



(48) 



Here Z and Z k are hermitian matrices, while S(0) has 
both real and imaginary parts and is a complex non- 
hermitian matrix. 

Next we compute the square root r k = yZk through 
the eigensystem of Z k . We thus have 

G k (w) = r k (w - r k (-/i + E(0) + e k - itu 2 B^y 1 r{^Q) 

We first solve the non-hermitian eigenvalue problem 



[r k (e k -fi + S(0))r k ] A k = AgE k 
At [r k (e k — fi + E(0))nj - ^k^, 

and compute the scattering rate in the eigenbase 

A£r k B k r k Ag = r k . 



(50) 
(51) 



(52) 



of the following functions 



Pi{z) = I dx 
P 2 (z,j) = / dx 



df\ 1 



(59) 



dx J x — z 

—) 1 ( 60 ) 

dx J \x — z + ix 2 "f\ 2 



If p = q, we have 



dx 



1 



dx J (x — z + ix 2 ^) 2 



Sp P — Q2(E kp ,T kp ) 



Rp P = P2(E kp ,T kp ) 



and if p ^ q we approximate 



Pi(E kp ) - P!(E kq ) 



I? 



E kp 
Pi(E kp )- 



- E kq 
Pi(E kq ) 



</;< 



E k p ~ E *q 



(61) 



(62) 
(63) 



(64) 
(65) 



Here we neglected the term proportional to x 2 T in the 
denominator, since the derivative of the fermi function 
constrains \x\ <C 1 and since the interband transition 
give subleading contribution to the Drude peak. 

A special care needs to be taken to compute the inte- 
grals Pi, P2 and Q2 to high enough precision and avoid 
divergencies. We give details on their evaluation in Ap- 
pendix [5] 



to get 



G k ( W ) 



r*A* 



Eu 



-Atr k 



(53) 



Here we used E koj = E k — iujT k . Next we insert Eq. (53) 



into ( 46 ) and we neglect the off-diagonal components of 
the scattering rate ( {T k ) pq ~ T kp 5 PA ), since the scatter- 
ing between quasiparticles is subleading at low tempera- 
ture. We thus obtain 



271-Fr 



^-ReJ] [Cp q S^ p - Dp q R 

k pq 



k 1 



(54) 



where 



C* = (A£r k v£r k AZ) qp (At<r k vZr k AZ) pq (55) 
Dp q = (Akr k v£r k A^) qp (A«\ k vZr k A«) pq (56) 
df\ 1 



c k = 

qp 



dx 



dx 



dx J (x - E kxp )(x - E kxq ) 
df\ 1 



J \txq ) 



(57) 
(58) 



The integrals and i? k have multiple poles and need to 
be treated by care. We first rewrite S k and i? k in terms 



VII. IMPURITY SOLVERS BASED ON 
HYBRIDIZATION EXPANSION 

The impurity solvers based on the hybridization ex- 
pansion have a long history and were often employed to 
solve the problem of a degenerate magnetic impurity in a 
metallic hoslP^^l In the past, most of calculations were 
limited to the lowest order self-consistent approximation, 
called the Non-crossing approximation (NCA). Recently, 
many generalization of the approach were studiecpSlSiJM] 
to overcome the difficulty of the NCA at low tempera- 
ture, below the Kondo temperature. It is well known 
that the NCA approximation fails to recover the Fermi 
liquid fixed point at low temperature and low energy. 
Typically there are three types of problems with NCA: i) 
the Kondo temperature is correct when only one type of 
charge fluctuations is dominant (like N — > N—l, which is 
equivalent to the limit of U — 00). When more than one 
charge fluctuation needs to be considered (N — >• TV + 1 
and TV — > N — 1) the Kondo temperature is severely un- 
derestimated and hence the Kondo peak is too narrow, 
ii) The asymmetry of the Kondo-Suhl resonance and its 
height is exaggerated in NCA. iii) At very low temper- 
atures T <C Tk an additional spurious peak at zero fre- 
quency appears. 

For DMFT applications, the problem iii) is not very 
severe, while the other two are. The first problem can 
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be corrected by a very moderate computational ex pense. 
Adding the first subleading Feynman diagr am d 25 * 54 !, 
named One crossing approximation (OCApEH cures the 
problem of the low energy scale. It also substantially im- 
proves the asymmetry of the Kondo peak as well as its 
width. Not surprisingly, in the context of DMFT, the 
OCA approximation gives correct critical U of the Mott 
transition in the Hubbard model, while NCA severely 
underestimates it. In c ontra st to other higher order con- 
serving approximations^^, the OCA approximation is 
relatively straighforward to generalized to the arbitrary 
impurity problem. Due its attractive features, OCA was 
used in many DMFT applications, such as unraveling the 
mixed valence state in PvP^. the coherence-incoherence 
crossover in Ce-115 materials^, the transport properties 
in titanides 1 -^, the a to 7 transition in Ce, etc. Com- 
pared to exact solution, as obtained by CTQMC, the 
OCA approximations typically gives very precise proba- 
bility for all atomic states^ (the histogram), quite pre- 
cise coherence scale, and the quasiparticle renormaliza- 
tion amplitude (the width of the Kondo peak), which is 
typically only slightly underestimated. At temperatures 
below the coherence scale, the OCA method, however, 
still suffers from slight overestimation of the height of the 
Kondo peak, and hence causality violation in the context 
of DMFT. Hence, the OCA approximation has to be used 
with care, especially in the systems with high coherence 
scale, and the systems with only moderate correlations. 

The OCA equations for the one band problem were 
given by many authors ^ 5 ! 54 ! , and their generalization to 
multiband situation was briefly discussed in the review 
Ref. [51 the generalized equations were however, not yet 
given, hence we will give them for the general multiorbital 
impurity problem, as relevant in the electronic structure 
calculations in section IVTIBl 

Recently, a renewed interest i n the hybridization ex- 
pansion arouse, once it was showrJ2i!22l fo^t the Feynman 
diagrams can be efficiently sampled by Monte Carlo im- 
portance sampling. The current implementation of this 
algorithm, as applied to realistic m ateria l problems, was 
discussed in plenty of detail recently22BU, and it will not 
be repeated here. 

Here we will rather outline an alternative Monte Carlo 
sampling approach, which was not yet discussed in the 
literature nor implemented. It is natural to ask if there 
exists an alternative regrouping of diagrams in Monte 
Carlo sampling, such that NCA approximation would be 
the lowest order contribution in the hybridization expan- 
sion, i.e., the two kinks approximation. We detail the 
method below in section IVII Al and show results of a 



simplified implementation, which truncates the sampling 
at a finite order (up to fifth order in hybridization). 



A. Towards Bold-CTQMC 



Ref. [Ml and Ref. EUJ). On the other hand, the OCA im- 
purity solver is very accurate in many correlated systems 
with narrow bands. For example, it gives correct criti- 
cal U in Hubbard model, correct Kondo scale in Kondo 
lattice model, etc. 

The current implementation of CTQMC is equivalent 
to pseudoparticle formulation of the expansion around 
the atomic limit, however, with bare pseudoparticle prop- 
agators. It is thus natural to expect that the dressed 
pseudoparticle propagators would make the algorithm 
more efficient, since the two kinks approximation is 
equivalent to NCA, and the four kinks approximation 
to OCA. 

The basic idea of the bold CTQMC algorithm is to 
sample the skeleton Feynman diagrams, with propaga- 
tors being dressed^. The Monte Carlo importance sam- 
pling samples all such diagrams, with the probability pro- 
portional to their Luttinger-Ward functional $. Hence 
contributions to all pseudoparticle self-energies can be 
straighforwardly sampled within this approach. 

Although the formalism of hybridization expansion on 
real axis was developed long ago (see for example ReflB^|) , 
its imaginary axis counterpart was not yet given. To 
our knowledge, the NCA equations have not yet been 
implemented on imaginary axis, because of the problems 
with diverging term in the projected Dyson equation (see 
Eq. ([79} below). 

In the hybridization expansion, the pseudoparticles are 
introduced to diagonalize the atomic part of the Hamil- 
tonian. The impurity problem is cast into the form 

H = ^2\m)E m (m\+'^2ekicl i Cki (66) 

m ki 
+ ^2 V kia \m)(m\fl\n)(n\c ki + h.c. 

mn,kai 

where we used completeness ^2 m \m)(m\ = 1 for atomic 
states I to). Each atomic state is represented by corre- 
sponding pseudoparticle ajj vacuum) = \m), and the 
completness of atomic basis gives a constraint for pseu- 
doparticles — > Ylm a m a m = Q = 1. The Hamiltonian is 
then given by 



* -ki 



(67) 



ki 



^2 V kla al n a n (m\fl\n)c ki + h.c. + X{Q - 1) 



mn,kai 

and the action is 



/d 
dTa) m (— + E m + X)a m 
in 

n'm' 

nn'mm' 

x / drdr / a^(r)a 7l (r)A QtJ g(T-r / )a^/(T / )o TO '(r / ) 



(68) 



The CTQMO^H solver is the most efficient exact 
solver for electronic structure problems (see for example 



where [F a *) 
XQ. 



Afi\ 



We also define H = H n 
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Any physical quantity has to be evaluated in the Q = 1 
subspace. This is achieved by letting A — > oo, to separate 
the spectra of Q = 0, Q = 1, Q = 2, • • • . Then we use 
the Abrikosov's trick to pick out the Q = 1 subspace. 
The expectation value, which we want to compute is 



(A) Q - 



Trr 



Tr Q=1 (e-^) ' 



(69) 



while accesible quantities are (A) — J2q Trq(^4e l3H )/Z. 
If operator A vanishes in the absence of impurity (in 
Q = subspace), the physical expectation value can be 
computed by 



(70) 



This is clear from expansion 



Z(A) = Tr Q=1 (^e 
Z(Q) = Tv Q=1 (. 



) + Tr Q=2 (Ae- 



Z = Tr Q=0 (e-^) + - 



i(2e 



-f3H -2p\ 



(71) 



Notice also that in the A 



limit 



(Q)eP x = 



Tr 



Q=i 



(72) 



can be used to obtain impurity free energy. 

In more general case, when (A) does not vanish in 
Q = subspace, Eq. (70) should be replace by (A)q = i = 
i im , (AQ1 

The Green's functions for pseudoparticles obey the 
Dyson equation, 



Gm, — 



1 



— A — En 



E m (w)' 



(73) 



where the energies of all pseudoparticles are shifted by 
A compared to atomic energies E m , due to XQ term in 
the Hamiltonian. In general, the Green's functions for 
pseudoparticles are off-diagonal. The states which corre- 
spond to the same superstate, defined in Ref. [531 obey 
a matrix analog of the above Dyson equation. However, 
here we will give equations for diagonal case, since the 
generalization is less transparent, but straighforward. 

The numeric limit of A — > oo is very untractable for 
computer. Since bold-CTQMC is implemented in imagi- 
nary time, we thus want to analytically project the pseu- 
doparticle equations on imaginary time axis. 

Before the limit A — > oo is taken, the pseudoparticlc 
Green's functions are given by 



G m (r) 



XT S~1 

On 

-XT 



(X) 
.(*) 



T > 
T < 



(74) 



The poles of the Green's function G m are at large fre- 
quencies, comparable to A, while G rn vanishes for x -C A. 



Hence G(r < 0) vanishes because f(x)G m (x) vanishes. 
We thus have 



G m (r <0)=0 
G m (r > 0) = e~ AT 



dx 



(75) 

r G' m (x + X). (76) 



This equations demonstrate the well known fact that the 
pseudoparticles can not propagate back in time. 

To derive a set of well posed projected equations, we 
introduce projected Green's functions, which remain well 
behaved in the limit A — > oo, and are used for numeric 
implementation 



G m (r) = e Xr G m (T) 



(77) 



Of course, these projected propagators vanish for r < 
0. The projected propagators are analogous to the well 
known projected functions on the real axis (see Ref. 
G m (x) = G' m (x + X)/f(-x) since 



G m (r) = 



r .f{-uj)G m {Lo) 



(78) 



is the usual r <-> lo transformation between the imaginary 
time and real frequency. 

Our goal is to write all equations in terms of projected 
G and analogous E functions, which do not contain A 
and are numerically well behaved. The problem however 
is that the projected quantities do not have fermionic nor 
bosonic character, and hence can not be represented on 
imaginary frequency axis. The Dyson equation Eq. ( 73 ) 
can be expressed in terms of projected functions by 



G(r) = T 



- (iuj — A)r 



- A - E - /,f dr'e^-^'t^ 1 ) 



(79) 



but its evaluation is far from straightforward. For conve- 
nience, we drop the index m from E m , G TO , and £ m . 

We need to evaluate this formula in the limit A — > oo. 
It is however not possible to perform the limit numer- 
ically because the exponential factors grow as exp(A/3) 
while the poles are in infinity on the real axis. 

For the implementation of the bold-CTQMC, it is cru- 
cial to find numerically tractable form of the projected 
Dyson equation. To this end, we perform expansion in 
powers of E, to get 



~ p -(kd-A)t 

g(t) = tJ2~. 



s 



iuj — X — E \ iu — X — E 

ILJ 

s 2 

+ {iuj-X-E) 2 



(80) 



where S — dT'e^ lu A ) r E(r'). The summation over 
imaginary frequency can now be performed, to obtain 
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G(r) 



E 



i d n 

n! dE" 



driE(ri) 



dT 2 S(r 2 ) 



dT„S(T„)e 



— E(t— n — r 2 r„) 



(81) 



Note that the limits of integration are constraint to 
the phase space of forward propagating pseudoparticles. 
Namely, the limit of A — s- oo does not allow the time 
difference in the exponent to be negative. 

To evaluate the projected Dyson equation in a stable 
way, we first evaluate the following moment-functions 



finite shift E.„ 



S„(t) 



1 

n! 



dT'E(r')e 



-E(t-t')( _ i\n 



(r-rT, (82) 



and then we convolve the moment-functions with £. 
Eq. (fSll) is hence implemented by 



The 



G(r) - -e-** + S 1 (T)-(i:*S 2 )(T) 

+ (£* (£*£ 3 ))(t) + (E*(E* 



where 



(S*Q)(r)= / dr'S(r - rOQCr') 
Jo 



(83) 
*5 n ))(r) 

(84) 



Note that all terms in the expansion have the same sign 
(note E < 0), hence the expansion converges quite fast, 
and we typically need between 30-50 terms for numeri- 
cally sufficient precison. 

Convolutions can be evaluated by standard method 
of Fourier transforms, or, they can be cast into the 
form of matrix multiplications, once the matrix E T]T / = 
E(t — t')cIt' is precomputed and used for all terms in the 
expansion. 

It is instructive to check the formula in two simple 
limits: i) £ = £ <5(t), evaluates to G(t) = - e -^+ s ") r ; 

ii) E = (To = const and E = evaluates to G(t) = 
-cosh(r^CTo"). 

The latter limit is very instructive because it shows 
that G can exponentially grow at low temperature and 
finite t. This is well known problem from implementing 
the NCA equations on real axis. To keep G(t) finite, 
and peaked around the origin on real axis (G(r) roughly 
constant in r), one needs to shift all pseudoparticle ener- 
gies E m — > E m + Xq to sufficiently positive energies, such 
that J2m~ G m (P ~ + ) = const, where const is of the 
order unity. Namely, in grand canonical ensemble, the 
pseudoparticle charge (Q) , defined in Eq. ( 72 ) , is 



(Q) = GmW - + ) = e-^ G '»^ ~ 0+ ) ( 85 ) 

m m 

indeed vanishes in the physical Q = 1 subspace. Once 
the projection is done, the physical quantities in Q = 1 
subspace are invariant with respect to shift of all pseu- 
doparticle energies by the same amount. If we introduce 



A 



E m + Aq (which is equivalent to 



A + Ao), charge (Q) will decrease for e while 
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will remain the same. 



the product (Q)e 
Similarly, all physical quantities are invariant, while the 
projected pseudoparticle quantities are not. Hence, for 
numerical stable evaluations, it is crucial to choose the 
shift Ao such that pseudoparticle propagators are finite. 
A large Ao will make them exponentially small, while van- 
ishing A will cause G m to diverge at j3. We thus need 
to fix the value of Ao properly. Two possible choices are 
Em G m(^ ~ 0+ ) = const or G m - gs ((3 - 0+) = const, 
where m — gs is the pseuodoparticle, which corresponds 
to the ground state of the atom. 

The basic idea for the bold-CTQMC is to sample 
self-energies for all pseudoparticles as well as the local 
Green's function. This is easiest to achive by defining 
the probabilty to be proportional to the absolute value 
of the Luttinger-Ward functional |$[G, A]|, and the self- 
energies then become 



G a fs — 



6$[G, A] 

1 <5$[G,A] 
W) SA fja 



(86) 
(87) 



where the first equation is contribution to the pseudopar- 
ticles self-energies, and the second is contribution to 
the real-electron Green's function (the impurity Green's 
function). 

The second identify might be less obvious, but it fol- 
lows from the fact that the impurity Green's function is 
the T-matrix for the conduction electrons 



1 



9k 



-E f 



: 9uSu,jk' + 9kiVki a G a pVk>ji3- 



ki,k' j 



We have seen above that G mm i carries a factor of e~ A/3 , 
and we will show below that $ also carries the same fac- 
tor e~ A ^, hence the pseudoparticle self-energy E mm ' is of 
the order of unity. On the other hand, the conduction 
electron self-energy E c is proportional to 5${G,A]/5A, 
and hence vanishes as e~^ x . Therefore both E c and G a p 
are propo rtio nal to e _/3A . The expansion of the equa- 
tion Eq. (88) in powers of e _/3A shows that i) conduc- 



tion electron propagator is unrenormalized in this the- 
ory (or equivalently the bare hybridization A appears 
in functional $[G, A]); ii) the impurity Green's func- 
tion, evaluated in the grand-canonical ensemble G a p is 
equal to (5$[G, A]/<5A, which vanishes as e~^ A . However, 
the physical quantities like the electron Green's function 
must be evaluated in Q = 1 subspace, using Eq. (70 1. 
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The resulting ratio is of order unity and is invariant with 
respect to shift of Ao, as explained above. 

The Luttinger-Ward functional $[G, A] for the lowest 
order contribution (two kinks), known under the name 
NCA, is given by 

$°[G,A] = / dTG mm ,{T)G n , n {(3-r)A a p{-T) 



x(F a ) nm (F^) m , n , (89) 

Note that if integration variable is shifted to r —> f3 — r, 
additional minus sign can appear. In case of regular 
fermions and bosons, this minus sign is automatically 
taken care of by the antiperiodicity of fermionic Green's 
functions G{{3 — r) = — G(— r). The pseudoparticle 
Green's functions however vanish at negative times, and 
one needs to add (3 to the negative argument, and add 
an overal minus sign when j3 is added to the fermionic 
Green's function. 

The corresponding pseudoparticle self-energies are 



>(r) = (-l) f 



6$°[G,A] 
5G n , n (f3-T 



(90) 



where (—1)-^ is +1 (-1) if n corresponds to pseudo-boson 
(pseudo-fcrmion) . Again, this minus sign is because neg- 
ative times are not allowed for pseudoparticlcs. 

Each pseudoparticle propagator carries an exponent 
e~ AAr , and the sum of exponents is always e~^ A . This 
holds for all diagrams composed of exactly one loop of 
pseudoparticles. These are the only diagrams that give 
contribution to the physical quantities. 



If we take out the exponential factors, the NCA func- 
tional takes the form 

<&° [G, A] = e~ px / drG mm ,(T)G n , n (P - t)A^(-t) 
Jo 

x(F") nm (F^) m ^,91) 
If we denote $[G, A] = e /3A $[G, A], we see that 
5$ <5$ 



S(r) 



E(r)e-^ 



SG(P-t) 5G{/3-t) 



hence 



6&[G, A] 

6G n , n (f3-r) 



(92) 



The projected $[G, A] has exactly the same form as 
$[G, A], we only need to replace G -> G. The NCA 
diagram hence becomes 



(93) 



$°[G,A] = / dTG mm ,(T)G n , n (/3-T)A afj (-T) 

m'n' 



From Eqs. (92) and rt93| it is clear that we achieved 



the goal of expressing all equations in terms of projected 
quantities, which do not depend on variable A, and are 
numerically well behaved. 

The projected second order diagram, which correspond 
to OCA approximation, is given by 



$ 1 {G,A]= / dr 4 / dr 3 / dr 2 / dnG mom ' o (n - r 4 + /3)G mim; (r 2 - T X )G m ^ (r 3 - r 2 )G TO3m ^ (r 4 - r 3 ) 
Jo Jo Jo Jo 

{F^) m , imo A a p{ Tl - r 3 )(F^) m>2 (F a '+) miroi A Q ^(r 2 )(^') m Gm3 (94) 
I 



The projected pseudoparticles vanish at negative times 
and are well behaved at positive times. For the purpose 
of properly evaluating the Feynman diagrams in time, 
we can extend them to negative times without any loss 
of generality. The pseudo-bosons hence become periodic, 
and the pseudo- fermions antiperiodic. The annoying mi- 
nus signs (—1)-^ can then be eliminated. However, the 
projected pseudoparticles can not be Fourier transformed 
to imaginary frequency, and they do not obey the usual 
Dyson equation, but rather a more complicated type of 
Dyson equations derived in Eq. (81). The pseudoparti- 



cles can be analytically continued to real frequencies, and 
all pseudoparticles satisfy fermionic-type of continuation, 
given in Eq. ( 78 ) . 



Finally, the Monte Carlo algorithm must generate any 



skeleton diagram of any order. The probability to accept 
the diagram is proportional to its |$[G, A]|. The contri- 
bution to pseudoparticle self-energy is then E mm /(r) = 
(sign($)/G m ' m (— t)), where () means the average in the 
Markov proces, where weights are proportional to |$|. 
Similarly, the impurity Green's function can be sam- 
pled by G a /j (t) = (sign(*)/A^ Q (-r))/(Q). The sam- 
pled self-energies will only be proportional to the exact 
self-energies. The renormalization factor can easily be 
found knowing the probability for NCA diagram, and its 
value. 

The requirement to sample the skeleton diagrams pro- 
hibits us to combine many diagrams into determinant 
of hybridization functions A, as it was achieved in the 
algorithm by Werner et.alpB Similar type of trick of 
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combining the diagrams into determinant of A's would 
substantially improve the efficiency of the algorithm, ft 
is however not clear how to eliminate non-skeleton dia- 
grams from determinant, and keep the updating formulas 
efficient. 

To test the above described algorithm, and to check 
its performance and convergence, we implemented a sim- 
plified version of the bold-ctqmc for the canonical An- 
derson impurity model. We sampled all diagrams up to 
certain order starting with first order (NCA), second or- 
der (OCA) and up to fifth order. The fifth order takes 
only minutes on a typical personal computer. We first 
found the topology of all diagrams of certain order, the 
prefactor and the sign of each diagram. In Fig. [2] we 
plotted the diagrams for the first few orders (second - 
$( 2 \ ... fouth - $( 4 )). We colored the dia grams accord- 
ing to their sign, positive with black and negative with 
red. There are four NCA diagrams, two OCA diagrams, 
8 third order diagrams (4 positive and 4 negative), 44 
forth order diagrams (24 positive and 20 negative), 320 
fifth order diagrams (128 positive and 192 negative). We 
evaluated exactly the NCA and OCA diagrams, and we 
used Metropolis algorithm to sample the time arguments 
for higher order diagrams. The probability for the ac- 
ceptance of a set of imaginary times was taken to be 
proportional to the value of the total |$( n ) (n, t 2 , ...t 2 „)|, 
hence at fifth order 320 diagrams were evaluated at each 
Monte Carlo step. While this algorithm can not be used 
at very high orders in perturbation theory due to expo- 
nential growth in the number of diagrams, its advantage 
is in large improvement of the sign problem. Namely, 
the diagrams of the same order and the same time argu- 
ments tend to cancel at higher orders. Since we evaluate 
all of them at each Monte Carlo step, the sign problem 
is almost completely eliminated. 

The non-interacting limit U — is the hardest case 
for the hybridization expansion algorithm, because the 
coherence temperature is infinite. Here we present test 
of the algorithm in the case of half-filled non-interacting 
Hubbard model on the Bethe lattice within DMFT. Wc 
want to emphasize that the algorithm becomes more ef- 
ficient and faster converging in strongly interacting limit 
U >> 0, a case which will be presented elsewhere. 

In Fig. [3](a) we show the impurity Green's function on 
imaginary axis (at l/T = 100) when the perturbation 
theory is truncated at certain order. We also display the 
exact result by the dashed line. While the NCA curve 
clearly deviates from the exact result, the higer order 
approximations are hardly distinguished from the exact 
curve on this plot. In Fig.[3]D we show separately the con- 
tributions to the Green's function from different orders in 
perturbation theory. As expected the contribution from 
the lowest two orders is large, while the higher order con- 
tributions are smaller. This shows why OCA approxima- 
tion is so successful in many realistic situations. The fifth 
order contribution is on average only 3 x 10 -3 , and never 
exceeds 6 x 10~ 3 . 

In Fig. [4] we zoom-in the exponential drop of the 




$(3)[G] = 





FIG. 2: All diagrams of the second, the third, and the 
fourth order in hybridization strength which contribute to the 
Luttinger-Ward functional. The pseudoparticle-propagators 
run across the ring, while the crossing lines stand for the 
hybridization A. The full line represents spin-up, and the 
dashed lines the spin-down hybridization. The black diagrams 
(both diagrams in $ (2) , first four in $ (3) , and first 24 in $ (4) ) 
give positive contribution to $, and the red give negative con- 
tribution. Some diagrams seems to appear multiple times. 
This is because different pseudoparticles appear in the ring. 
Since we do not use different line for each pseudoparticle, some 
diagrams seem equivalent. However, it is very straighforward 
to deduce the pseudoparticle propagators knowing the type 
and the direction of the conduction electron propagators. 



Green's function at short times. We see that the conver- 
gence with the perturbation order is very encouraging. 

For efficiency of the bold-ctqmc, it is important to 
monitor the sign of each individual diagram. In Fig. [5] we 
show separately the contribution to the impurity Green's 
function from the diagrams with positive <f> and those 
with negative <&, together with the sum of the two. At 
the third order, the sum is around 70% of the positive 
contribution, while at the forth and fifth order, the sign 
drops to 0.2 and 0.07, respectively. As explained above, 
the current implementation of the method, which groups 
together all diagrams of a certain order in perturbation 
theory, does not have a substantial minus sign problem. 
However, this method becomes expensive at high orders, 
and thus one needs to resort to sampling of individual di- 
agrams, which can be performed to arbitrary high order. 
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FIG. 3: (a) The comparison of the finite order perturbation 
theory result with the exact impurity Green's function, (b) 
The contributions to the impurity Green's function up to the 
fifth order, plotted separately order by order. 



FIG. 5: The three panels show the contribution to the impu- 
rity Green's function at 3rd, 4th and 5th order in perturbation 
theory. We show separately the contribution from the terms 
with positive $ and the terms with negative $. 




FIG. 4: The same as in Fig. [3] but we zoom in the short time 
behaviour. 



In the latter case, there will be a minus sign problem, as 
estimated here. 



B. The One crossing approximation 

In this section we will give the most general formulas 
for the One crossing approximation, and we will explain 
the crucial steps in implementing the algorithm. 

We start with lowest order approximation, which is the 
Non-crossing approximation. When evaluating these di- 
agrams, we have to consider only two Hilbert subspaces 
of constant N at once, i.e., N and N + 1. The first 
step is to compute all eigenvalues and eigenvectors of the 
atom in the subspace N and N + 1. We then group to- 
gether the atomic eigenstates, which are degenerate, i.e., 



have the same atomic energy E m . In the next step we 
check which of these degeneracy's survive in the presence 
of the crystal field environment (impurity hybridization 
A), and which off-diagonal propagators need to be con- 
sidered. We evaluate the following matrix elements 



b 



/ £deg ,(a ,a' )£deg and A Q 



(F a ') b2f (F^) fbl (95) 



>¥=o, 



Here b runs in the Hilbert subspace of N and / in 
the Hilbert supspace of N + 1. The matrix elements 
(F a ) bf = (b\f a \f) and (F"t) /6 = (/|/t| 6 ) where f a is 
electron destruction operator. The sum runs only over 
the / states which are degenerate and over one electron 
states a which are also degenerate and for which A aa > 
is nonzero in the considered crystal field symmetry. The 
resulting matrix elements Cb 2 ,b t have the same symmetry 
as the propagators of the pseudoparticles Gb 2 bi- Clearly, 
in high symmetry crystal environment, most of the off- 
diagional matrix elements vanish and the degeneracy of 
Gbb is high, but in low symmetry environment and in the 
broken symmetry state, many of the off-diagonal propa- 
gators become crucial. 

Once the symmtry of the propagators is known, we 
determine all nonvanishing bubbles (NCA diagrams) and 
the matrix elements for each bubble. The NCA matrix 
elements are 



°Oi6 2 /l/2 



E 



{F a ') b2h (F^) hbl m 



(/l>/2),(bl,b 2 ),(a,"')Sdeg 



where we sum only over degenerate states /, b and degen- 
erate crystal field components a. The Luttinger-Ward 
functional and the self-energy corrections are depicted in 
Fig. 6j We associate a factor (F a1< ) fi to each vertex that 
marks the creation of electron in bath a. Accordingly, we 
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FIG. 6: The NCA Luttinger-Ward functional and the self- 
energies within NCA. 
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FIG. 7: The Luttinger-Ward functional for the One Crossing 
Approximation (OCA). 



add a factor (F a )y for each vertex of electron anhilation. 

In the next step, we precompute the matrix elements of 
the one-crossing diagrams, which are depicted in Fig. [7] 
Here we need to select three different Hilbert subspaces: 
N - 1, N, and N + 1 to compute 

jjaa' 130' 

flf2hf4bib 2 a 1 a 2 

^2(F^) b2h (F a ') fiai (F^) a2h (F^) f2bl (97) 

deg 

Here b, f, a run over the states with N —1, N and N + 1 
number of particles, respectively. We add only the most 
important crossing corrections, for which the particle 
number N is in the Hilbert subspace of the ground state 
of the atom. We also select fi to be only the ground state 
multiplet of the atom, or the atomic states with energy 
very close to the ground state energy. We compute the 
matrix elements C and D only once in the DMFT self- 
consistent loop and we save them into the input file for 
OCA impurity solver. The matrix elements C, D do need 
to be updated in the outer LDA+DMFT charge loop. We 
typically update them every three to four charge steps, 
since the relative crystal field splittings usually change 
very little during LDA+DMFT iterations. The atomic 
energies E m change much more (due to the chemical po- 
tential shift), and need to be updated at every step. 

The NCA diagrams on the real axis can be evaluated 
with conventional techniques, and after the projection, 
they take the following form 



5W")= E -( Fa ')»2h(F^) hbl J 'ff(y)K a ,(y)G flh (u; + y) 



(98) 



(99) 



bi b<20LOi' 



]T J dye-^(F a ') b2fl (F^) hbl Gl b2 (y)Gl f2 (y + ^ 



bib 2 fif: 



where G — ImG. The pseudoparticle propagators G and 
the pseudoparticle self-energies are related by the Dyson 
equation. The Eq. 67 shows that G = 1/ '{omega — E — 
A-S). 

Many of the pseudoparticle propagators and hybridiza- 
tion functions are degenerate, hence in practice we do not 
need to sum over all possible b, f and a indices, but we 
rather use the precomputed matrix elements Cj^j^, 
which make sure that no equivalent diagram (a diagram 
which has the same frequency dependence) is not com- 



puted multiple times. 

To take care of the diverging exponential factors, we 
work with the projected quantities G(uj) = G (cj) //(— u>) 
and S(w) = S w), as explained above. The 

pseudoparticles have typically very sharp almost diverg- 
ing structure near the treshold energy, which is not easy 
to Fourier transform. Hence we can not use the Fourier 
transform for convolutions. We rather cast the above 
equation into the form for matrix multiplication, for 
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which fast linear algebra packages such as BLAS, exist. 
We use the logarithmic mesh to resolve the fine structure 
of the pseudoparticle green's functions. 

It is important to realize that the number of baths a 
is quite small (of the order of 2(2L + 1) for correlated 
orbital of angular momentum L), while the number of 
atomic states is much bigger. Hence we precompute the 
integral and the first moment of functions A (w)f{uS) 
and of A (uj)f(—uj) for all aa' . Within trapezoid rule, 
the values and the first moments of these quantities are 
enough to compute the above convolutions with matrix 
multiplications on any given mesh. 

To see that, lets consider an arbitrary convolution 



(100) 



C(z) = / g(x)f(x — z)dx 



Here the function g{x) is defined on a certain mesh {xi}, 
on which it is well resolved, i.e., g(xi) = g^. The function 
f(y) is defined on another mesh {yj, i.e., /(?/,) = /j. 
The convolution can be safely calculated on the union of 
both mashes {xi,yj + z}. One of the meshes should be 
shifted for z, thus for each outside frequency, a different 
union of the two meshes should be formed and only then 
the convolution can be safely evaluated. This is very time 
consuming and not done in practice. 

When a certain / function needs to be convolved with 
many other functions (like A (co)f(u>) in our example 
above), we use the followin trick. We first precompute 
the integral and the first moment of the function 



We then calculate the convolution without building a new 
inside mesh. Let's use the mesh {x;} which resolves func- 
tion g. Then, in the spirit of trapezoid rule, we can lin- 
early interpolate g between the points 



C(z) = 



9i+i 



-[x — Xi) f{x — z)dx. 

(104) 

This integral can be expressed by the above defined func- 
tions. To show that, let us rewrite the convolution and 
expressed it by the new function (/),; which is defined on 
the same mesh as g and with which the covolution is a 
simple scalar product 



C{z)=Y,9i 



x i + 1 -z 



Xr l + l — Z — U 



Ci — z Xi 
~ Z Z + U — Xi-i 



f(u)du 



f{u)du 

Xi _ 1 — Z Xi Xi—i 

= J2gi{fhzdhi. (105) 



fate) = 



f(u)du 
uf(u)du 



(101) 

(102) 
(103) 



Thus (f) lz is 



(f)i 



'(xj+i - z)[Fi(x i+1 - z) - Fi(xj - z)] - F 2 (x l+ i - z) + F 2 (xj - z) 
{x i+ i - Xi-i)(x i+ i - Xi) 
(x i _ 1 - z)[Fi(xj - z) - F^x^i - z)} - F 2 (x l - z) + F 2 (x i _ 1 - z)\ 
(x l+ i - Xi~i)(xi - Xi-i) 



(106) 
(107) 



Hence the convolution of / with many functions g m can 
be computed at once C(m, z) — g rn * / by the following 
matrix product C(m,z) = Y^i9mi{f)izdhi. 

Once the NCA contributions are evaluated, we add 



the second order diagrams, which correspond to OCA 
approximation and are depicted in Fig. [7J They take the 
explicit form 
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£ b2bl (u;) = - 



fihf 3 f4.aia 2 ce0a'f}' 



(108) 



/ d ^.f(y) A 'w(y)Gf 3 fA^ + y) J ^f{x)A n aa ,{x)G hb {w + x)G aia2 ( w + x + y ) 



£ 020 » = - E (^)feA(i fa ')/ 1 .i(^)^(^)w 1 x 

fif 2 f 3 f4bib 2 af3a'0' 



(109) 



| ^/(-l/)A^(l/)G /8/4 (a;-y) | ^/(-^^(^^(a; - :r)G tl6 > - a; - y) 



E [(^)w3(f a ')M 1 (^ , ) a2 / i r t )w 1 + (F a ) b2fl (F ')f 2ai (F a \ 2 f 3 (F^)u bl 



f 3 f4.aia 2 bib 2 af)a'fl' 



J d ^f(-y)Al a ,(y)G blb2 (uj - y) J ^ f{x)A^(x)G aia2 (w + x)G hft (u + x - y) 



(110) 



47( £ ) = - E 

act' fif 2 f 3 f4bib 2 aia 2 



{F l3 ') b2 f 3 (F a ')f 4ai (F^) a2 f 1 (F a ^)f 2bl + (F a ') b2fl (F ')f 2ai (F a ^) a2h (F^) f4bl ](idl) 



1 f f dx a 

s" A (Q)/(-e) J ^ ^ J V /(x)AaQ,(x)Im{G6lh2(y)G/l/2(x + y)}Im{G/3/4(e + 2/)Gaia2(£ + a; + y)} 



In practice, we do not sum over all /, b and a indices. 
As explained above, we precompute the matrix elements 

F > f^f 2 f 3 f 4 b 1 b 2 a 1 a 2 f° r tnc mos t important processes. We 
take only the low lying atomic states into account (only 
/'s which are part of the ground state multiplet or with 
energy very close to the ground state). We also take 
into account the degeneracy of all atomic states and the 
degeneracy of baths a, in order to avoid computing the 
equivalent diagram mutliple times. Finally the convolu- 
tions for the OCA approximation can also be cast into 
the form of matrix multiplication, once the first moment 
and integrals of a few functions are precomputed. 



VIII. THE ANALYTIC CONTINUATION 
METHOD 

The Monte Carlo impurity solvers are implemented on 
imaginary axis where the quantity being sampled is real 
and many times even "sign-free" . The results obtained in 
this way are exact, except for the statistical noise. How- 
ever, even a tiny statistical error on imaginary axis pre- 
cludes the analytic continuation by Pade type of meth- 
ods. The standard method, to overcome the difficulty of 
the singularity of the kernel, is the Maximum Entropy 
Method (MEM). The basic idea of this method is to 
find a function on the real axis, which is very close to 
Monte Carlo data on imaginary axis (within statistical 
error), and is smooth function on real axis, locally not 



very different from a chosen model function. This ap- 
proach works very well for analytical continuation of the 
Green's function G(r) to obtain spectral function on real 
axis, i.e., to solve the integral equation 



0(t) 



//(-,) 



e TX A(x)dx 



for A(x). 

Knowing the spectral function, it is however not possi- 
ble to obtain the momentum resolved spectra, or optical 
conductivity, or transport coefficients. To compute these 
properties, it is essential to analytically continue the self- 
energy, rather then the Green's function. The self-energy 
of correlated materials is however very hard to analyti- 
cally continue with maximum entropy method, because 
the self-energy typically has very sharp feature or even 
poles, which separate the low energy part of the spec- 
tra (the quasiparticle peak) from the high energy part 
of the spectra (the Hubbard bands). Due to the max- 
imum entropy method requirements of smoothness, the 
analytically continued self-energy at low energy is typ- 
ically polluted with the near-by poles, which appear in 
the self-energy at the intermediate energy. 

A successful analytic continuation method for self- 
enery needs to met the following conditions: 

• imaginary axis self-energy is equal to Monte Carlo 
data within the statistical error 

• real axis self-energy function must be locally 
smooth 
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• the power-expansion around zero frequency should 
match the quantum Monte Carlo data on both, real 
and imaginar axis. 

While the first two conditions are met by MEM, the last 
is not. 

We developed an alternative method, which mets the 
above conditions and was very successfully used in com- 
bination with CTQMC for pnictidesP, cupratesP, V0 2 
and other materials. Although the method has many pa- 
rameters, which needs to be choosen appropriately, we 
can always check its accuracy by recomputing the spec- 
tral function of the lattice, using analytically continued 
self-energy, and comparing the spectral function to the 
maximum entropy continued spectra. 

We expand the self-energy in terms of modified Gaus- 
sians £, and we add a polynomial function around zero 
frequency 



S M (z) 



E 



The modified Gaussians 
C"(E n ,u)= 1 



c n £(E n ,z) + f (z) 



,-b 2 /<L-(lag(u/E n )/bf 



b\E n \^7T 



(112) 



(113) 



have a unique peoperty that they are peaked at E„, with 
the width of approximately E n , while they exponentially 
vanish at zero frequency. They are asymmetric with slow 
decay away from zero and very fast decay towards zero 
frequency. We choose the modified Gaussians centered 
on a logarithmic mesh of E n — ±7rTw" with w ~ 1.5. 
The modified Gaussians functions were used in connec- 
tion with constructing the NRG spectral function^. We 
typically take the parameter b to be ~ 0.8. 



Since the modified Gaussians all vanish at zero fre- 
quency, we add a polynomial function around zero fre- 
quency. The coefficients of the polinomial are determined 
by fitting the imaginary axis self-energy, i.e., 

S(w„) = S + (-61 + iai)w n + (-02 - ih)^ (114) 

which can be analytically continued to 

E(w) = S + (<Zi + ibi)uj + (a 2 + ib 2 )u 2 (H5) 

The polynomial has to drop-off sufficiently fast at high 
frequency, hence we choose the following function 



fo M = 



So + ubi + uo 2 b 2 ) I (1 + (uj 2 b 2 /T 2 ) 2 ) 



S T 2 /(u 2 + T 2 ) 



where the upper choice is made for metals and the lower 
choice for insulators and very bad metals. In the FL 
regime, we have b\ <SC 1, |Sq| oc Z 2 ir 2 T 2 . The coefficient 
T is determined by the condition f (oj = 1) <C 1 

For speed, we precompute C(E n , iui) and C (E n ,co) by 



C{E n , Z ) = - l ~ ( dxC {E - X) . (117) 

7T J Z — X , 

Similarly, we also precompute fo(iu>) and /q(uj). Also the 
integral of the functions I n = J dxC (E n ,x) and I = 
J dxf (x) are precomputed. 



The coefficients c„ in expansion Eq. (112 1 are deter 



mined by minimizing the following functional 



X = E I S m(^„) - S QMC (iu;„)| 2 + a 1 |/(S M ) - I D \ 2 + a 2 |S M (0) - S | 5 

uj n £sampled 



-a 3 



rf s M (Q) 

du) 



(01 + ibx) 



«4 



rf 2 S M (0) 
duj 2 



2{a 2 +ib 2 ) 



(118) 
(119) 



Here uj n in the first term runs over the imaginary fre- 
quencies which are sampled by QMC (and not over the 
analytically added tail). The second terms imposes the 
correct value of the integral of the self-energy. The inte- 
gral of the expansion ( |112[ ) is 

I(Em) = c n I n + 7 , 

n 

which needs to match the l/(oj n ) tail of the QMC data 
I D =ir lim uj n E QMC (uj n ) 



Finally, the last three terms ensure that the value, and 
the first two derivatives of the analytically continued self- 
energy at zero frequency match the derivatives on imag- 
inary axis. For minimization, we use the L-BFGS-B al- 
rorithm of Ref . 66 . 

IX. CERIUM a-7 TRANSITION 

To test our implementation of DFT+DMFT within 
Wien2K method, we show in Figs. [8] and [9] results for 
cerium a to 7 transition. 
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FIG. 8: Total and partial density of states of elemental cerium 
metal in both phases, a and 7 phase. We used OCA impurity 
solver. 



FIG. 10: The hybridization function of the j z = 5/2 subshell 
within LDA and within DMFT in both phases. 
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FIG. 9: The same data as in Fig. [8] but obtained by continu- 
ous time quantum Monte Carlo solver, and analytical contin- 
uation method. 



At a temperature less than 600 K and pressure less 
than 20 kbar, elemental cerium undergoes a transition 
between two isostructural phases: a high pressure phase 
or a phase and a low pressure 7 phase. In a-Ce the 
/ electron is delocalized while in 7-Ce the / electron is 
localized. The transition is well accounted for by phe- 
nomenological Kondo Volume Collapse picture™™. 

We treat only the Ce 4 — / electrons as strongly corre- 
lated thus requiring full energy resolution, while all other 
electrons such as Ce spd are assumed to be well described 
by the GGA. We choose U = 5.5 eV and J = 0.68 eV for 
the Coulomb interaction. The value of U was obtained 
by constraint DFT calculation^ and J was computed us- 
ing the atomic physics program of Ref. and reduced 




FIG. 11: Optical conductivity of a-Ce and 7-Ce within 
LDA(Wien2K)+DMFT method. Note the shoulder in a-Ce 
conductivity, which is due to excitations across the two quasi- 
particle peaks (4/ : 5/2 and 4/ : 7/2) clearly visible in figure 
[8j and also measured by experiment of Ref. \ST\ 



by 30% to account for the screening in the solid. Both 
phases of Ce have fee unit cell with quite different vol- 
umes, V a = 28.06A 3 and V 7 = 34.37A 3 . The results were 
converged with 5000-k points, we use the GGA functional 
for the DFT part and use OCA and CTQMC impurity 
solver to solve the auxliliary impurity problem. 

The results in Fig. M (obtained by OCA) and Fig. [9] 
(obtained by CTQMC) are practically identical and very 
similar to previous LDA(LMTO) +DMFT results^. One 
can clearly see the broad quasiparticle peak in a-Ce, split 
by the spin-orbit coupling ~ 0.3 eV. The lower peak has 
mostly 5/2 character and the upper peak mostly 7/2- 
character. The system is in good Fermi liquid regime 
at the temperature of 150 K used in the calculation. 
The second phase with larger volume is in local moment 
regime with no visible Kondo peak at the Fermi level, 
but enhanced Hubbard bands. 
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It is instructive to examine the hybridization function 
A = lu — Ei mp — £ — 1/G as computed by LDA and self- 
consistent DMFT (see Fig.[lO]). It turns out that in Ce, 
the low energy hybridization function is substantially re- 
duced compared to its LDA value. The two large peaks 
at — 0.4eV and 0.7 eV are absent in DMFT hybridization. 
Since the coherence scale is exponential function of hy- 
bridization, the coherence scale is lower in DMFT than it 
would be in so called one-shot DMFT. It is known from 
the early days of the Kondo volume collapse theory^, 
that the LDA hybridization in a one-shot calculation was 
too big and had to be renormalized by phenomenological 
parametei^. DMFT reduces the hybridization through 
the collective screening effects and hence is able to give 
correct coherence scale of the problem. 

Finally, let us show optical conductivity, as imple- 
mented in LDA ( Wien2K) +DMFT method. The overal 
agrement with previous LDA+DMFT resultfP is very 
good. The new computational results are in even slightly 
better agreement with experiment of Ref. |HZ1 since they 
both clearly display a shoulder around 0.3 eV in a-Ce, 
which we can now clearly identify as excitations across 
the split quasiparticle peak. The splitting is due to spin- 
orbit coupling in Ce. 



X. HEAVY FERMION 115 MATERIALS 

The heavy fermion 115 materials have a chemical for- 
mula CeAIn 5 , where X is either Co, Rh or Ir. They crys- 
tallize in layered tetragonal structure shown in Fig. 12 
composed of Ce-In layers and A-In layers. 

At high temperature, the low energy electronic states 
are composed of mainly the broad spd bands of In and 
Ce. The Ce-4/ electrons are localized and their spectra 
is mostly contained in Hubbard bands, which are more 
than 2eV away from the Fermi level. These electrons 
behave as local magnetic moments. As the temperature 
is reduced, the moments combine with the conduction 
electrons to form a fluid of very heavy quasiparticles, 
with masses that are two or three orders of magnitude 
larger then the mass of the electrons. 

The low temperature physics of 115 materials is very 
puzzling. The heavy fermion physics comes primarily 
from the Ce-In layer. Indeed, the related material Celn 3 
has only the Ce-In layers (no X-ln layer), and also dis- 
plays a similar heavy fermion properties with supercon- 
ductivity at very low temperature. However, 115 mate- 
rials are very sensitive to the substitution of the transi- 
tion metal ion in the X-In layer although Co, Rh and 
Ir ions have the same valence (they are isovalent). In- 
deed the three 115 materials have dramatically different 
low energy properties: CeCoIn 5 is a superconductor with 
T c ~ 2.3 K, CeRhIn5 is antiferromagnet with T/v ~ 3.5 K, 
while Celrlns is superconductor with T c of only 0.4 K. A 
fundamental question arises: Why are the low energy 
properties of 115 materials so different? 

A hint to the resolution of this problem was given in 




FIG. 12: Crystal structure of CeXIns. Red, yellow and gray 
spheres correspond to Ce, X , and In atoms, respectively. 



Ref. EH1 where the DFT+DMFT calculation for Celrln 5 
indicated that the Ce 4/ electrons hybridize stronger to 
the out of plane In-p electrons, than the in-plane In p 
electrons. Here we carried out the DFT+DMFT calcula- 
tion for all three 115 materials and we show the difference 
in electronic structure between the three materials. We 
used the code based on LDA-LMTO code of Ref. [3U as 
well as the new LAPW code based on Wien2Kpil code. 
The results obtained by our DFT+DMFT method in the 
two codes are almost indistinguishable. For the impu- 
rity solver, we used both OCA (described above) and 
CTQMCP. The analytic continuation of CTQMC re- 
sults was performed with the method described in chap- 
ter Em 



Fig. 13 K shows the total density of states (DOS) and 



the partial Ce-4/ DOS for all three materials at low tem- 
perature of 7K. The transition metal ion DOS is peaked 
around binding energy 2eV, where the difference of DOS 
is large. The partial Ce-4/ DOS of the three compound 
is very similar, except at the very low energy. Fig. |13} 3 
zooms-in the low energy part of the spectra. We see that 
Celrlns compound has the largest quasiparticle peak, the 
CeCoIns follows, while the CeRhlns has substantially 
smaller quasiparticle peak at the same temperature of 
7K. 

Our view on the localization- itinerancy in 115 materi- 
als is sketched in Fig. 14 Rh compound is most local- 



ized, while Ir compound is most itinerant. Co compound 
is similar to Ir compound, but slightly less itinerant than 

Ir-115. 

It is well known from the pressure experiments^^ 
that Rh compound is more localized then Co compound. 
Namely, under pressure of 1 GPa the Rh compound 
becomes superconducting, and at pressure of ~ 2 GPa 
reaches similar maximum T^ as is the maximum Tc 
of Co compouncPSJ. Hence the pressure of the order of 
GPa sufficiently increases the Ce-4/ hybridization that 
it overcome the difference between localization of the 
electrons in the two compounds. Experimentally it is 
a bit less clear what is the relation between Ir and 
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FIG. 13: Total density of states (full lines) and partial Ce- 
4/ density of states (dashed lines) for CeCoIns, CeRhlns and 
Celrlns materials. The lower pannel show the low energy part 
of the Ce-4/ density of states for all three compounds. We 
used OCA solver. 



Co compound, since both compounds are superconduc- 
tors at low T. Ir compound has somewhat smaller spe- 
cific heat coefficient in normal state than Co compound 
(mo lK 2 ) for Ir-115 versus 1000 mJ/(mo/K 2 ) for 
Ir compound has also somewhat lower resis- 
tivity in the normal stated. Moreover, nuclear quadrupol 



(750 m i I (mo 
Co-lisp^. 



resonance (NQR) measurements of 1/(T\T)°°- suggest 
that Ir-compound might be more itinerant than other Ce- 
compound. Indeed pressurizing the CelrlnJ^II along the 
crystallographic c-direction, which increases itinerancy^!, 
decreases Tc- Furthermore, it was shown that Cd-doping 
acts as reverse pressure in 115's®. Since higher Cd- 
doping is necessary for appearance of antiferromagnetic 
phase in Celrlns than in CeCoIns, this is also suggestive 
of more itinerant nature of Ir-compound. 

Our results are thus consistent with the resistivity 
experiments^, NQR experiments^! and recent pressure 
experiments^, and indicate that Ir compound is on the 
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"good" Fermi liquid superconducting magnetically ordered 
Ce-ln(l) 6.246 6.164 6.222 

Ce-ln(2) 6.183 6.202 6.194 

angle 45-0.59 45+0.35 45-0.26 



FIG. 14: The sketch of the itinerancy/localization of the 
three 115 compounds. In our view, the Ir compound is most 
itinerant, while the Rh compound is most localized. The Co 
compound is not localized enough to develop magnetic order 
at low temperature, while it is nor a good metal. It is thus 
conceivable that it would show tendency towards supercon- 
ductivity. This phenomena is however beyond our current 
theoretical method - the single site DMFT calculation. We 
also show the bond distances between Ce and In atoms, and 
the angle between Ce and out-of-plane In atom. None of 
these parameters can explain the actual order of the com- 
pounds, hence the structure itself can not explain the trend 
of localized to itinerant transition in these compounds. 



itinerant side of the phase diagram. Hence the low super- 
conducting transition temperature might be connected 
with too itinerant nature of carriers. 

We further analize the difference in itinerancy by plot- 
ting the hybridization function at zero frequency A(o> = 
0)ll resolved in crystal field basis. The 14 dimensional 
matrix of hybridizations has a 6 dimensional j = 5/2 
component and a 9 dimensional j = 7/2 component. 
The itinerancy (the quasiparticle peak) is almost entirely 
from the j — 5/2 component, hence we will not analize 
j = 7/2 part. The degeneracy of the 5/2 shell is lifted in 
tetragonal crystal environment and hybridization splits 



into r v 



F+ 
1 7 



and Tq components. 



The Tq corresponds 
to j z = ±1/2, while and correspond approxi- 
mately to j z = 3/2 and j z = 5/2, respectively. In Fig. 15 
we plot the hybridization (— Im[A(0)]/7r) in polar coor- 
dinates with Ce atom in the center and In atoms around. 
The plot is the cut in xz direction. The three dimensional 
orbitals that correspond to the three crystal fields are 



plotted in the lower pannel of Fig. 15 together with the 
real space positions of In and Ce atoms. In this decom- 
position, hybridizations and are pointing towards 
out of plane In (I112) and in-plane In (Ini), respectively. 
The third component, Tq is pointing towards transition 
metal ion. 

When comparing hybridization of the three 115 com- 
pounds, the Ir compound has all three components of 
the hybridization larger than the other two compounds. 
In Co-115 all three hybridizations are slightly smaller, 
while in Rh-115 all three hybridizations are substantially 
smaller. 

Furthermore, comparing the strength of the three com- 
ponents of the hybridization, one can notice that in Ir 
compound the Ty component, pointing towards out-of 
plane In, is largest. This is consistent with the exper- 
imental finding of Oeschler et alW^ that the Griineisen 
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FIG. 15: The Ce-4/ Weiss field hybridization function \A"\/tt 
decomposed into crystal field components of tetragonal field. 
All quantities are in units of eV. The upper plot shows the 
2D projection of the three relevant orbitals, while the lower 
pannel shows their 3D shapes (blue dots mark the position of 
the in-plane In atoms, while the red dots the position of the 
out-of-plane In atoms). The radial extend of the orbitals in 
the polar plot of the upper panel is proportional to the value 
A /n at zero frequency. The full/dashed/dotted lines corre- 
spond to Celrlns/CeCoIns/CeRhlns. While all three compo- 
nents of hybridization F^T, F^", and Tg, are largest (smallest) 
in Celrlns (CeRhlns) compound, F^" takes the largest value 
and also changes more than the other two components. 



parameters in c-direction is 2.5 times bigger that in a 
direction, resulting in larger effective coherence scale in 
c-direction. 

In Co-compund the and Tq components have simi- 
lar strength, while Ty is smaller, hence the hybridization 
in c-direction is still more important than in afc-plane, 
consistent with Griineisen parameter measurements^!. It 
was shown in Ref. EH that the double peak structure of 
the optical conductivity is directly related to the strength 
of the two hybridizations. The hybridization gap in one 
part of the momentum space is larger, and is primarily 
due to out of plane In, and the hybridization gap in some 
other part of momentum space, controled mainly by the 
in-plane In, is smaller, resulting in double peak structure 
of the mid-infrared optics peak. Optical measurements 
on CeCoIn5 of Singley et demonstrated very clearly 
that the mid-infrared peak is split into two peaks, one at 
250cm -1 and one at 630cm _1 , which can hint towards 
substantial difference in the two types of hybridization. 



Finally, in contrast to Ir and Co compound, Rh com- 
pound has largest Tq hybridization, followed by and 
IV. 



0.05 



J. 0.03 
0.02 
0.05 
n 0.04 



Jt 0.03 

< 

0.02 



— Co 

— Rh 

— Ir 













-0.4 -0.2 0.0 0.2 0.4 



FIG. 16: The frequency dependence of the three most im- 
portant hybridization functions A (u})/n in all three 115 com- 
pounds: CeCoIns, CeRhlns and, Celrlns. The frequency is 
in units of eV. 

In Fig. [16] we show the frequency dependent hybridiza- 
tion function — ImA(w)/7r to demonstrate that the retar- 
dation effects in heavy fermion materials are very non- 
trivial and that the buildup of the quasiparticle peak in 
spectral function usually results in a sharp peak in hy- 
bridization, on the background of the depleted region of 
hybridization. The peak is sometimes called the collec- 
tive hybridization, because it arises from the lattice ef- 
fects. Namely, the Ce-4/ electrons on neighboring atoms 
also become delocalized, enhancing the hybridization at 
low energy. However, the spd-electrons need to screen 
many Ce-4/ moments, and therefore the effective spd hy- 
bridization is actually slightly reduced, resulting in deple- 
tion away from the Fermi level, sometimes called Kondo 
hole. 

Our results demonstrate that the degree of itinerancy 
is controlled by the collective hybridization, encoded into 
the Weiss mean field hybridization A(u>) within DMFT. 
But what is the origin of the difference between the three 
compounds? In Fig. [14] we show the parameters of the 
lattice structure, namely the Ce-In(l) distance, the Ce- 
In(2) distance and the angle between the Celn3 plane 
and out of plane In (R12). From these numbers, it is 
clear that none of the three quantities follows the trend 
of itinerancy. Hence the difference in the lattice structure 
is likely not the key element. 

To demonstrate that the difference in the lattice struc- 
ture is not the driving force, we performed the DMFT 
calculation for the three compounds using the same lat- 
tice structure of Celrln 5 . The results were very similar 
to the results plotted in Fig. [l3j with only slight increase 
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in itinerancy of Rh compound. This demonstrates that 
the chemistry of the transition metal ion (difference be- 
tween 3d, Ad and 5d orbitals) is the driving force of the 
itinerancy, and not the diffrence in the crystal structure. 
The latter are the secondary effects. 




FIG. 17: Temperature dependence of the magnetic moment 
of the comensurate AFM Neel state in CeRhlns. 

Since CcRhlns remains in local moment regime down 
to very low tempeture of the order of the RKKY inter- 
action, it is worth trying to stabilize a magnetic solution 
within DMFT. To this end, we doubled the unit cell and 
allowed the comensurate antiferromagnetic ordering with 
the wave vector (1/2,1/2,1/2). Experimentally, the or- 
der is a helical spiral with wave vector (1/2,1/2,0.298) 
and Tc of 3.8 K. The broken symmetry solution can be 
stabilized below T ~ 3 K as shown in Fig. [ITJ The mag- 
netization has a typical mean field form, as expected for 
a theory with spatial mean-field character like DMFT. 

An interesting question is how does the large moment 
antiferromagnetic solution change the emerging quasipar- 
ticle peak. We have shown in Fig. 13 that even in more 
localized CeRhIn 5 a peak starts to develop at the Fermi 
level by decreasing temperature, hence coherence starts 
to develop at quite high temperature similar to the other 
two compounds. However, the height of the quasipar- 
ticle peak is smaller and the scattering rate of Ce-4/ 
orbital (imaginary part of the self-energy) is higher in 
CeRhIn 5 . The long range order state develops from a 
state with a partially screened moment. In Fig. [18] we 
show the density of states of the two phases, the para- 
magnetic state and the Neel state. The latter has no 
quasiparticle peak left and only a very broad background 
of the / spectral weigh remains at the Fermi level. The 
lower panel of Fig. [18] compares a very coherent quasi- 
particle peak of Celrlns with the partially screened state 
of CeRhIn 5 above T^ ee i and in the ordered state below 
Tpf ee i, to emphasize the dramatic difference in the den- 
sity of state at low energy. Because the full coherence of 
quasiparticles is not reached to very low temperature in 
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FIG. 18: Total and partial Ce-4/ density of states for Rh- 
115 below and above the AFM transition. Above the Neel 
temperature, there is a signature of Kondo effect, wich par- 
tially screenes magnetic moment at elevated temperatures, 
even though the system develops the long range order below 
3K. The quasiparticle peak is however much smaller than the 
same peak in Celrlns material. Once in the ordered state, the 
quasiparticle peak dissapears. 



CeRhlns, and the non-local RKKY interaction is strong 
enough, it interrupts the formation of coherent quasi- 
particles. Within DMFT, this is reflected in two stable 
solutions of DMFT equations, the paramagnetic and the 
magnetic solution. We note that we did not prove the 
stability of the magnetic solution compared to the para- 
magnetic solution, because this would require a compar- 
ison between free energies, a task beyond our current 
capabilities. However, our experience from model calcu- 
lations suggests that when the magnetic DMFT solution 
can be stabilized, it usually has lower free energy than 
the nonmagnetic solution. 
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XI. CONCLUSION 

In the first part of the article, we discussed in de- 
tail the implementation of DFT+DMFT in full potential 
methods. We defined the central object of the DMFT, 
the local Green's function using a projection operator. 
We showed that the projector used in LDA+U imple- 
mentations leads to non-causal DMFT equations and 
that the straightforward projection to the solution of the 
Schrodinger equation within the Muffin Tin sphere leads 
to spectral weight loss. We suggested an alternative pro- 
jection that resolves these shortcomings. 

We sketched the algorithmic steps within an implemcn- 
tion of DFT+DMFT in the full potential methods, using 
a formulation which avoids the ambiguities of downfold- 
ing or Wannier orbital construction. Hence, the kinetic 
energy operator and electron density are not approx- 
imated by a tight-binding parameterization, which al- 
lowed us to carry out a charge density self-consistent cal- 
culation. 

In the second part of the article, we concentrated on 
impurity solvers based on the hybridization expansion. 
We derived the equations for the bold continuous time 
quantum Monte Carlo (CTQMC) method, which samples 
the dressed propagators, as opposed the bare propagators 
sampled in current CTQMC methods. We showed a few 
test results for simplified implementation of the method. 
In this part of the article we also gave detailed formulas 
for the impurity solver called the One-crossing approxi- 
mation, which can be viewed as the four kink approxi- 
mation within the bold CTQMC. 

Finally we give details on a new analytic continua- 
tion method, which can continue the self-energy from the 
imaginary to the real axis. This step is crucial when com- 
puting the response functions within DMFT, as done in 
section [Vl| for transport coefficients. 

In the third part of the article, we presented the test 
results of our DFT+DMFT implementation on a classical 
problem of strong correlations, the isostructural transi- 
tion of elemental cerium from its 7 phase at high tem- 
perature to its a phase at low temperature. 

In the last part of the article, we applied the 
DFT+DMFT method to a group of heavy fermion com- 
pounds, namely Celrln 5 , CeCoIn 5 and CeRhIn 5 , collec- 
tively dubbed the Ce-115s. Although the isovalent sub- 
stitution of a transition metal ion does not substantially 
alter the Ce-In planes, which are believed to be respon- 
sible for the heavy mass in these compounds, the ground 
state properties of these materials are very different. 

We analyzed the electronic structure of the three Ce- 
115 materials and showed that the Ce-4/ electrons in 
CeRhlns are more localized that those in the other two 
115 compounds, in agreement with experiments. Be- 
low 3K, an antiferromagnetic DFT+DMFT solution in 
CeRhlns is stable, while CeCoIns and Celrlns remain 
paramagnetic (the AFM solution is not stable) down to 
the lowest temperature T = 1.5 K explored in our calcu- 
lation. 



The hybridization in Celrlns is very anisotropic with 
the largest component pointing towards the out-of-plane 
In. The hybridization is slightly smaller in CeCoIn 5 , 
hence we believe Celrlns to be more itinerant than the 
other two compounds. 

We speculate that the reason CeCoIns exhibits the 
highest superconducting Tq is due to the fact that it is 
at the border between itinerancy and localization, while 
Celrlns is on the itinerant side of the phase diagram 
and CeRhIn5 is on the localized side. The position of 
CeRhIn 5 in the phase diagram is clear from the pres- 
sure experiments, while the position of Celrln5 is less 
obvious. We believe that recent uniaxial pressure exper- 
iments^ confirm our view, since the c-axis compression, 
which makes Celrln 5 more itinerant, decreases the super- 
conducting T c . 
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Note added: As the writing of this work was be- 
ing completed, we became aware of a related work of 
M. Aichhorn et al. (arXiv: 0906.3735) also reporting on 
an implementation of LDA+DMFT in a LAPW code. 
However, in contrast to our implementation, the authors 
used downfolding method to obtain a tight-binding model 
Hamiltonian, and hence could not compute electronic 
charge self-consistently. In the process of downfolding, 
the authors used projection P 1 defined in section [n] 
Furthermore, the impurity solver did not take into ac- 
count the full Coulomb interaction Eq. ( 28 1 . Only the 



density-density part of the interaction was considered by 
M. Aichhorn et al. (only the z component of the Hund's 
coupling) which allows substantial simplification of the 
impurity solver, but leads to improper description of the 
multiplet structure of the correlated atoms. 



Appendix A: Complex Tetrahdron Method 

The formulas for tetrahedron integral in case of com- 
plex eigenvalues are very similar to the case of real eigen- 
values. However, a special attention needs to be payed 
to choose the right branch-cut in logarithms, such that 
all terms in the sum are causal. 

First step in tetrahedron method consists of dividing 
the first Brillouin zone into tetrahedra which fill up whole 
space. Each thrahedra has four corners. The energy is 
thus interpolated e = e\ +a(e 2 — £1) + &( £ 3 + c(£4 — 
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£i), where a, b and c run between and 1 when visiting 
corners of tetrahedra. 

For the Green's function we need integral of the form 



E 



w(k,^)C k 



(Al) 



where 



ilv(x,2/) = -j/ 2 { w 4 [log(x) -log(x + y)] +log(x + y) 



+ u 3 + -u 2 - u] 



(A7) 



and for the electron density and the chemical potential 
we need 



E 



did- 



Cu 



w - £k 



^m(k,w)C* k (A2) 



The integral is first written as the sum over all tetra- 
hedra and the integral in the interior of tetrahedra: 

k t K t fc 4 =l 

The latter is evaluated analytically using linear interpo- 
lation inside the volume of the tetrahedra for both the 
nominator and denominator 

/>1 pX—c rl—b—c 

w(ki,u) = 6 dc db dax (A4) 

Jo Jo Jo 

(1 - a - b - c)5 ku i + a8k t ,2 + b5 ki ,3 + c8k u A 

u)-£x- a(e 2 - ei) - &O3 - £i) - c(e 4 - £1) 

Here we used a short notation = 

The integrals are analytic and a closed expression for 
computing the green's function is 



UJ — £,- 



lv (w - £,-,£j - £*) 



where 



lv(x, y) = ^ |l - I [log(x + y) - log(x)] | (A5) 

and I 7^ i,j means I 7^ i and I ^ j. Notice that only 
log(x + y) and log(x) can appear in lv(x, y) (not log(y)) 
to ensure causality. Namely, imaginary part of all Si is 
strictly negative, hence the expression lv (lj — Ej, Ej — £$) 
contains log(<J — Ej) and log(w — £j), which both have 
imaginary part in the interval [0,7r]. 

Similarly, the formulas for the integral over frequency 
JiT] 2 w(ki,uj)du} are 



wi(ki, u>2, oji) 



ilv(w 2 - £j,ej - Si) 



Zlilli Zfij (A6) 
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and u = x/y 

Appendix B: Transport integrals 

To compute the transport coefficients, we need to eval- 
uate to high precision the following integrals 



df\ 1 



Pi (z) = J dx I - 



/ x — z 
df\ J 



Q2(2,7) = J dx ^- 



dx J |x — z + ix 2 "/^ 
df\ 1 



dx J (x — z + ix 2r y)' 



(Bl) 
(B2) 
(B3) 



The integrals need to be carefully implemented and 
special care needs to be taken for the two case: a) \z\ 1 
and b) \z"\ < 1 and \ j\ < 1. 



The first integral of Eq. (Bl I is computed numerically, 



except in the following cases 

J -(l/z + co/z 3 + ci/z 5 + c 2 /z 7 + c 3 /z 9 ) \z\ > 10 
n(Z) \w (z') +l n^(z') |/|«1 



where 



t — X 



is precomputed on a fine mesh and interpolated using 
cubic spline interpolation. The constants Cj are 



7vr 4 



317T 6 
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Co = y Cl = 15"' C2 = ^r 



C3 = 
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(B4) 



The second integral of Eq. (B2) is computed numeri 



cally, except in the following cases i) \z | <C 1, 7I -C 1: In 
this limit it becomes ^2(^,7) ~ J^\d^( z ); ") M > 14: 
In this case, the power expansion in terms of \z\ 2 is per- 
formed and all terms are analytically evaluated. 

Similarly we treat integral Eq. (B3|. For \z | <C 1, 
| 7 | « 1 we approximate Q 2 (z,i) ~ T ^*{§l) 

and for \z\ > 8 we perform the power expansion in terms 
of z 2 and analytically evaluated the resulting integrals. 
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